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ABSTRACT 

We introduce a dust model for cosmological simulations implemented in the moving- 
mesh code arepo and present a suite of cosmological hydro dynamical zoom-in simula¬ 
tions to study dust formation within galactic haloes. Our model accounts for the stellar 
production of dust, accretion of gas-phase metals onto existing grains, destruction of 
dust through local supernova activity, and dust driven by winds from star-forming 
regions. We find that accurate stellar and active galactic nuclei feedback is needed to 
reproduce the observed dust-metallicity relation and that dust growth largely dom¬ 
inates dust destruction. Our simulations predict a dust content of the interstellar 
medium which is consistent with observed scaling relations at 2 = 0, including scal¬ 
ings between dust-to-gas ratio and metallicity, dust mass and gas mass, dust-to-gas 
ratio and stellar mass, and dust-to-stellar mass ratio and gas fraction. We find that 
roughly two-thirds of dust at z = 0 originated from Type II supernovae, with the 
contribution from asymptotic giant branch stars below 20 per cent for z > 5. While 
our suite of Milky Way-sized galaxies forms dust in good agreement with a number of 
key observables, it predicts a high dust-to-metal ratio in the circumgalactic medium, 
which motivates a more realistic treatment of thermal sputtering of grains and dust 
cooling channels. 
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1 INTRODUCTION 

Dust in the interstellar medium (ISM) exists alongside gas- 
phase metals and alters the dynamic and spectroscopic prop¬ 
erties of galaxies (Calzetti, Kinney & Storchi-Bergmann 
1994; Silva et al. 1998; Dey et al. 1999; Calzetti et al. 2000; 
Netzer et al. 2007; Spoon et al. 2007; Melbourne et al. 2012). 
The surfaces of dust grains play host to a range of chemi¬ 
cal reactions that subsequently influence the behaviour of 
the ISM and impact star formation (Hollenbach & Salpeter 
1971; Mathis 1990; Li & Draine 2001; Draine 2003). Ad¬ 
ditionally, observations suggest that dust is a significant 
contributor of metals in the circumgalactic medium (CGM; 
Bouche et al. 2007; Peeples et al. 2014; Peek, Menard & Cor- 
rales 2015). Understanding the life cycle of dust grains, in¬ 
cluding their production in asymptotic giant branch (AGB) 
stars and supernovae (SNe; Gehrz 1989; Todini & Ferrara 
2001; Nozawa et al. 2003; Ferrarotti & Gail 2006; Nozawa 
et al. 2007; Zhukovska, Gail & Trieloff 2008; Nanni et al. 
2013; Schneider et al. 2014), growth via accretion of gas 
particles in the ISM and coagulation with other dust par¬ 
ticles (Draine 1990; Dominik & Tielens 1997; Dwek 1998; 
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Hirashita & Kuo 2011), destruction via thermal sputtering, 
collisions with other dust grains, and SN shocks (Draine & 
Salpeter 1979a,b; McKee 1989; Jones, Tielens & Hollenbach 
1996; Bianchi & Ferrara 2005; Yamasawa et al. 2011), and 
other physical processes, is important in accurately mod¬ 
elling dust evolution. 

Even at high redshift, galaxies can form substantial 
masses of dust. Far-infrared and submillimeter observations 
show that some dust-rich radio galaxies and quasars out 
to z ~ 7 have dust masses greater than 10 7 M 0 (Fan et al. 
2003; Bertoldi et al. 2003; Hughes, Dunlop & Rawlings 1997; 
Venemans et al. 2012; Casey, Narayanan & Cooray 2014; 
Riechers et al. 2014; Watson et al. 2015). Two notable ex¬ 
amples are SDSS J1148+5251, a z — 6.4 quasar with an in¬ 
ferred dust mass of (3.4^J'J®) x 10 8 M 0 (Valiante et al. 2009, 
2011), and A1689-zDl, a z = 7.5 ± 0.2 galaxy with a dust 
mass of 4 x 10 7 M 0 and a dust-to-gas ratio comparable to 
the Milky Way value (Watson et al. 2015). To produce such 
large dust masses at high redshift, average SNe must yield 
roughly 1 M 0 of dust, a quantity larger than the amount of 
dust SNe have been observed to condense (Todini & Ferrara 
2001; Sugerman et al. 2006; Dwek, Galliano &, Jones 2007; 
Lau et al. 2015). However, these dusty examples may not 
be representative of typical high-redshift galaxies. There is 
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recent evidence of actively star-forming galaxies at z > 6.5 
with little dust obscuration (Walter et al. 2012; Kanekar 
et al. 2013; Ouchi et al. 2013), and observations suggest 
that dust extinction is generally decreased in low luminos¬ 
ity and high redshift systems (Bouwens et al. 2012). Our 
current understanding of dust formation at high redshift is 
therefore still quite incomplete. 

Significant variation in dust properties also exists at 
low redshift. The Galactic dust-to-gas ratio is roughly 10“ 2 
(Gilmore, Wyse & Kuijken 1989; Sodroski et al. 1997; Zubko, 
Dwek & Arendt 2004) and several times larger than the 
values observed for the Large and Small Magellanic Clouds 
(Pei 1992; Gordon et al. 2014; Roman-Duval et al. 2014; 
Tchernyshyov et al. 2015). In contrast, the metal-poor local 
dwarf galaxy I Zwicky 18 has been estimated to have a dust 
mass of no more than 1800 M© and a corresponding dust- 
to-gas ratio in the range of (3.2 — 13) X 10 —6 , several orders 
of magnitude below the typical values expected for larger 
systems (Fisher et al. 2014). Additionally, observations and 
models of the ISM indicate that the nature of dust can differ 
among individual chemical species (Wilms, Allen & McCray 
2000; Kimura, Mann & Jessberger 2003; Jenkins 2009). 

The current sample of galaxies with reliable dust esti¬ 
mates has grown in recent years, driven by programs like 
SINGS (Kennicutt et al. 2003; Draine et al. 2007), the Her- 
schel Reference Survey (Boselli et al. 2010), the Herschel 
ATLAS (Eales et al. 2010), KINGFISH (Kennicutt et al. 
2011), and the Dwarf Galaxy Survey (Madden et al. 2013). 
Various trends and scaling relations involving dust and host 
galaxy properties have been observed from this data. There 
is a positive correlation between dust-to-gas ratio and metal- 
licity (Vladilo 1998; Draine et al. 2007; Galametz et al. 2011; 
Remy-Ruyer et al. 2014; Zahid et al. 2014). However, the de¬ 
tailed behaviour of the dust-metallicity relation is unclear. 
Observations indicate that the dust-to-gas ratio is reduced 
at low metallicity, possibly due to effects in the interstellar 
radiation field that limit dust growth and enhance destruc¬ 
tion processes (Remy-Ruyer et al. 2014). Even over small 
metallicity ranges, there is pronounced scatter in the dust- 
to-gas ratio. 

The observed dust-metallicity relation implies a dust- 
to-metal ratio that is fairly constant across a range of galaxy 
morphologies and histories. However, the nature of the dust- 
to-metal ratio at high redshift and low metallicity is uncer¬ 
tain. Recent work using gamma ray bursts has yielded dust- 
to-met al ratios fairly consistent with the Local Group, even 
in low metallicity systems (Zafar V Watson 2013; Sparre 
et al. 2014). This would require SNe to be efficient produc¬ 
ers of dust or grains to grow rapidly in the ISM (Mattsson 
et al. 2014). Separate analysis of gamma ray burst damped 
Lyman-alpha absorbers suggests a non-universal dust-to- 
metal ratio, with low metallicity environments producing 
dust less efficiently than spiral galaxies (De Cia et al. 2013). 
Even in the Milky Way, observations of strong gas-phase 
depletion (Roche V Aitken 1985; Sembach & Savage 1996; 
Jenkins 2004, 2009) contrast the expectation that dust de¬ 
struction outpaces stellar injection of dust (Barlow 1978; 
Draine & Salpeter 1979a; McKee 1989; Dwek & Scalo 1980; 
Jones et al. 1994; Jones, Tielens & Hollenbach 1996). Un¬ 
derstanding the balance between gas-phase metals and dust 
in Milky Way-sized systems requires more work. 

A number of other empirical scaling relations have 


emerged, tying observed dust masses to various galactic 
properties. These include relations between dust-to-stellar 
mass ratio and gas fraction (Cortese et al. 2012), dust-to- 
stellar mass ratio and redshift (Dunne et al. 2011), dust ex¬ 
tinction and stellar mass (Zahid et al. 2014), dust mass and 
gas mass (Corbelli et al. 2012), dust mass and star formation 
rate (da Cunha et al. 2010), and dust surface density and 
radial distance (Menard et al. 2010; Pappalardo et al. 2012). 
Observational data have also yielded initial estimates of the 
dust mass function for low and high redshift (Dunne et al. 
2000, 2011). These scaling relations provide constraints on 
galaxy formation models that include a treatment of dust. 

A variety of numerical models have been used in pre¬ 
vious work to better understand how dust evolves in a 
galaxy. These include one- and two-zone models (Dwek 1998; 
Lisenfeld & Ferrara 1998; Hirashita & Ferrara 2002; Inoue 
2003; Morgan & Edmunds 2003; Calura, Pipino & Mat- 
teucci 2008; Valiante et al. 2009; Gall, Andersen & Hjorth 
2011a; Yamasawa et al. 2011; Asano et al. 2013a; Zhukovska 
2014; Feldmann 2015), semi-analytic methods (Somerville 
et al. 2012; Mancini et al. 2015), and more recently first 
smoothed-particle hydrodynamical simulations resolving lo¬ 
cal dust variations (Bekki 2013, 2015a). These models in¬ 
clude processes like the formation of dust during stellar evo¬ 
lution, dust growth and destruction in the ISM, radiation 
field effects, and dust-enhanced molecular formation. Many 
of these one-zone models have addressed the formation of 
dusty high-redshift quasars and have indicated that if dust 
is unable to condense efficiently in SNe, a high star forma¬ 
tion or accretion efficiency is needed to produce large dust 
masses at z > 5 (Morgan & Edmunds 2003; Michalowski 
et al. 2010; Gall, Andersen & Hjorth 2011b; Gall, Hjorth & 
Andersen 2011; Michalowski 2015). 

Numerical dust models have also focused on the evolu¬ 
tion of the dust-to-gas and dust-to-metal ratios, often used 
in observations to characterise galaxies. The dust-metallicity 
relation is reproduced in many one-zone models (Issa, Ma- 
cLaren & Wolfendale 1990; Lisenfeld & Ferrara 1998). These 
models predict a dust-to-metal ratio that is independent of 
metallicity in galaxies where SNe are the primary producers 
of dust (Morgan & Edmunds 2003) and whose radial gradi¬ 
ent in galactic discs can be used to estimate interstellar dust 
growth (Mattsson, Andersen & Munkhammar 2012). The 
evolution of the dust-to-metal ratio may vary significantly 
with time, with the dust-to-metal ratio at z > 2 possibly 
just 20 per cent of the present value (Inoue 2003). Interest¬ 
ingly, galaxies are predicted to switch from low to high dust- 
to-metal ratio when crossing a critical metallicity thresh¬ 
old that enables efficient ISM dust growth (Zhukovska, Gail 
& Trieloff 2008; Inoue 2011; Asano et al. 2013a; Feldmann 
2015). Previous one-zone dust models applied to the for¬ 
mation of a single galaxy show present day dust-to-metal 
ratios of roughly 0.5 (Dwek 1998) and 0.9 (Calura, Pipino 
& Matteucci 2008), with results significantly dependent on 
the strength of dust growth and destruction mechanisms in 
the ISM. These models also find depletion roughly constant 
across all chemical species. However, the biggest weakness of 
these one-zone models is their lack of spatial resolution, lim¬ 
iting their ability to make predictions about the distribution 
of dust within a galaxy. 

Cosmological simulations provide a better means to un¬ 
derstand how dusty systems can form at high redshift and 


© 2016 RAS, MNRAS 000 , 1-28 


Dust Formation in Milky Way-like Galaxies 3 


how their dust content, both in the ISM and CGM, evolves in 
time. Additionally, simulations of full cosmological volumes 
can provide the sample size of galaxies needed to corrobo¬ 
rate observed scaling relations involving dust. Recent cos¬ 
mological hydrodynamical simulations have suggested that 
heavily dust-attenuated galaxies can form even at z ~ 7, 
with ultraviolet and optical colors consistent with a model 
having a dust-to-metal ratio of 0.4 (Kimm & Cen 2013). 
Motivated by the observed reddening of background quasars 
by foreground galaxies in Sloan Digital Sky Survey (SDSS) 
data (Menard et al. 2010), smoothed-particle hydrodynam¬ 
ical simulations found that half of this reddening signal is 
attributable to dust more than 100 h -1 kpc from the closest 
massive galaxy (Zu et al. 2011). Radial gradients of the dust- 
to-gas ratio in galactic discs appear to be steeper for larger 
galaxies (Bekki 2015a). Even in simulations where dust is 
not directly treated, radiative transfer post-processing can 
be used to infer dust extinction (Jonsson 2006; Narayanan 
et al. 2010; Hayward et al. 2013; Yajima et al. 2014). The 
dynamics of grain sputtering in SN shocks can also be stud¬ 
ied using tracer particles (Silvia, Smith h Shull 2010). While 
some studies have investigated the impact of feedback mech¬ 
anisms, such as winds driven by SNe, on the cosmological 
evolution of dust (Zu et al. 2011; Bekki 2015a), much work 
still remains. 

In this work, we incorporate essential dust physics into 
a moving-mesh simulation code and, for the first time in 
studies of dust, use a large sample of zoom-in cosmological 
initial conditions. Compared to previous dust models us¬ 
ing one-zone methods and idealised initial conditions, our 
approach has the ability to resolve the structure of dust 
within a galaxy and the impact that local feedback processes 
have on dust evolution. While some smoothed-particle hy¬ 
drodynamical simulations have treated dust without using 
zoom-in cosmological initial conditions, our highest resolu¬ 
tion cosmological run offers improved mass resolution. 

This paper is structured as follows. In Section 2, we dis¬ 
cuss the theory and numerical implementation of our dust 
model in the context of a larger galaxy formation frame¬ 
work. Details of the cosmological initial conditions for our 
simulations are provided in Section 3. We analyse the impor¬ 
tance of physical feedback processes, individual dust model 
components, and variations in simulation initial conditions 
in Section 4. Finally, in Section 5, we summarise our re¬ 
sults and discuss the implications of our findings on future 
observations. 


2 METHODS 

We employ cosmological simulations to track the evolution 
of dust and its effect on galaxy dynamics. We briefly re¬ 
view the core components of the galaxy formation model 
currently implemented in the moving-mesh code AREPO 
(Springel 2010), which have been detailed in prior work (Vo- 
gelsberger et al. 2013) and used for various cosmological 
studies (Vogelsberger et al. 2012; Genel et al. 2014; Torrey 
et al. 2014; Vogelsberger et al. 2014a, b), and then describe 
the new dust physics we have added. AREPO uses a dynamic 
Voronoi tessellation to solve the equations of ideal hydro¬ 
dynamics with a finite-volume method. A second-order Go¬ 
dunov scheme is used in conjunction with an exact Riemann 


solver to compute fluxes between cells. Additionally, gravita¬ 
tional and collisionless physics have been implemented using 
methods similar to the TreePM scheme in the smoothed- 
particle hydrodynamics code GADGET (Springel 2005). 

2.1 Galaxy Formation and Feedback Mechanisms 

The current galaxy formation model in AREPO accounts 
for a number of physical processes, including primordial 
and metal-line gas cooling, stellar evolution and subsequent 
chemical enrichment of the ISM, black hole formation, and 
stellar and active galactic nuclei (AGN) feedback, which to¬ 
gether yield a galaxy stellar mass function in good agree¬ 
ment with observations over 0 < z < 3 (Torrey et al. 2014; 
Genel et al. 2014). The star formation prescription stochas¬ 
tically creates star particles in dense regions of ISM gas, 
with masses distributed according to a Chabrier (2003) ini¬ 
tial mass function (IMF). As stars evolve off of the main se¬ 
quence, mass is recycled to the neighboring ISM. This chem¬ 
ical enrichment routine follows nine elements (H, He, C, N, 
O, Ne, Mg, Si, and Fe). We adopt elemental mass yields 
and recycling fractions for AGB stars from Karakas (2010), 
SNe II from Portinari, Chiosi & Bressan (1998), and SNe la 
from Thielemann et al. (2003). We employ a stellar lifetime 
function from Portinari, Chiosi & Bressan (1998). Galactic 
outflows are modelled using wind particles launched from 
star-forming ISM gas and temporarily decoupled from hy¬ 
drodynamics to emulate a galactic fountain driven by stellar 
feedback (Springel h Hernquist 2003). We include the mi¬ 
nor feedback modifications detailed in Marinacci, Pakmor V 
Springel (2014), which alter radio-mode AGN feedback and 
adopt warm galactic winds. 


2.2 Dust Evolution 

A wide variety of dust models have been used in recent 
galaxy formation simulations, though essentially all mod¬ 
els track both the formation of dust during stellar evolution 
and its subsequent evolution in the ISM. The most simplis¬ 
tic of these models take a one-zone approach, solving a set 
of coupled ordinary differential equations for the time evo¬ 
lution of the total mass of gas, metals, and dust within a 
galaxy (Hirashita & Ferrara 2002; Inoue 2003; Morgan V 
Edmunds 2003; Calura, Pipino h Matteucci 2008; Valiante 
et al. 2009; Gall, Andersen V Hjorth 2011a; Yamasawa et al. 
2011; Asano et al. 2013a; Zhukovska 2014; Feldmann 2015). 
Other work assumes a two-component galaxy model, con¬ 
sisting of a bulge and disc region and allowing for local 
dust properties to be studied as a function of radial dis¬ 
tance from the galactic center (Dwek 1998). More recently, 
Bekki (2013, 2015a) has performed smoothed-particle hy¬ 
drodynamical simulations that treat dust locally and offer 
improved spatial resolution. 

In this work, we focus on the dominant dust production 
and evolution mechanisms. We leave additional processes, 
such as the catalysis of molecular hydrogen on dust grains 
and the effect of interstellar radiation fields, to future efforts. 
We track the mass of dust in each chemical species within 
individual gas cells. Dust is passively advected between gas 
cells when solving the hydrodynamic equations in each time 
step, in essence adopting a two-fluid approach with dust 
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fully coupled to gas. However, dust is not implemented as a 
strictly passive tracer: as shown below, dust impacts metal¬ 
line cooling and in turn star formation and gas dynamics. 
Alternative treatments, including “live” dust particles de¬ 
coupled from gas and subject to frictional or radiative forces 
(Kwok 1975; Draine & Salpeter 1979b; Barsella et al. 1989; 
Theis & Orlova 2004; Bekki 2015b), may be pursued in fu¬ 
ture work. 

2.2.1 Dust Production via Stellar Evolution 

During the stellar evolution process described in Section 2.1, 
a certain amount of the mass A Mi of species i evolved by 
stars and returned to neighboring gas cells in some time step 
is assumed to condense into dust, with the exception of H, 
He, N, and Ne. The remaining metal mass exists in the gas 
phase. In the framework below, we follow the approach used 
by Dwek (1998) and other subsequent works that track in¬ 
dividual chemical elements. Most notably, we adopt a differ¬ 
ent functional form for the amount of dust produced during 
mass return from AGB stars than from SNe. Additionally, 
we make a distinction between AGB stars with C/O > 1 in 
their stellar envelope, which are expected to produce car¬ 
bonaceous solids (e.g. graphite or amorphous carbon), and 
those with C/O < 1, which are thought to form primarily 
silicate dust (Draine 1990). In essence, the equations below 
describe elemental mass yields for dust, computed as a func¬ 
tion of the yields for overall metals returned by stars. 

For AGB stars with a carbon-to-oxygen ratio of C/O > 
1 in their returned mass, the amount of species i dust pro¬ 
duced is given by 


and (3) is slightly reduced from that given in Dwek (1998). 
We have found that this minor modification is necessary for 
oxygen, whose dust production is tied to the ejecta of heav¬ 
ier elements, in order to avoid the formation of more oxygen 
dust than total oxygen mass returned. We use different stel¬ 
lar nucleosynthesis yields than in Dwek (1998), and even 
though there is considerable uncertainty in the condensa¬ 
tion efficiencies, we demonstrate in Section 4.3 that our re¬ 
sults are insensitive to moderate changes in these prefactors. 
The range of condensation efficiencies that we explore is de¬ 
scribed in Section 3.2. 

Within each gas cell, we separately track the mass of 
dust produced by AGB stars, SNe la, and SNe II, motivated 
by a desire to understand the dominant channels of dust pro¬ 
duction at high redshift. Additionally, we refer to carbon 
grains as graphite dust, and the dust in remaining species 
as silicate dust. While a simplification that fails to distin¬ 
guish the full diversity of compounds that comprise dust 
grains, including polycyclic aromatic hydrocarbons, this di¬ 
vision has been adopted in similar studies and is a step to¬ 
wards analysing the underlying ISM chemistry (Tielens & 
Allamandola 1987; Dwek 1998; Weingartner & Draine 2001). 

2.2.2 Interstellar Dust Growth 

The mass of dust in the ISM may increase over time, as gas- 
phase elements collide with existing grains (Draine 1990). 
To model the accretion of dust grains in the ISM, we follow 
the prescription of Dwek (1998) and Hirashita (1999) and 
in every time step compute each cell’s instantaneous dust 
growth rate 


AA/^diist — 


^GB’C/o > x (AMc - 0.75 AMo) 

0 


if i = C 
else, 

(i) 


where ^ GB,C ^° > 1 is the carbon condensation efficiency for 
AGB stars with C/O > 1, discussed below in more detail. 
Similarly, for AGB mass return with C/O < 1, the mass of 
species i dust formed is 


AM* dust — 


E 


rAGB,C/0 < 1 


j=Mg, Si, Fe 
.AGB,C/O < 1 AM 


A Mj / fij 


if z = C 
if i = O 

else, 

( 2 ) 


where m is the mass in amu and y^ GB ’ c /° < 1 j s the conden¬ 
sation efficiency for species i in AGB stars with C/O < 1. 
Finally, the mass of dust for element i produced via SNe 
ejecta is 


{ Sg N AM c ifz = C 

10 £ <5f N A Mj/fij if * = O (3) 

j=Mg, Si, Fe 

<5f N AMi else, 


where the dust condensation efficiency of element i for SNe is 
Sf N , which may differ from the corresponding efficiency for 
AGB stars. Equation (3) and its condensation efficiencies 
are used for both SNe la and II, consistent with all afore¬ 
mentioned models. Additionally, we note that the numerical 
prefactor for the calculation of AMo, dust in Equations (2) 


I dMj ; dust \ _ f dust \ f -Mj,dust 

V ) g \ -Mumetal / \ Tg 

where Mi, dust is the cell’s mass of species i dust, summed 
over components originating from AGB stars, SNe la, and 
SNe II, Mumetal is the corresponding species i metal mass, 
and r g is a characteristic growth timescale. Note that the 
first factor in parentheses induces a growth rate that de¬ 
pends on the local dust-to-metal ratio and slows the accre¬ 
tion rate as gas-phase metals are condensed into dust (Dwek 
& Scalo 1980; McKee 1989). 

Previous studies have explored accretion timescales de¬ 
pendent upon local gas density and temperature (Yozin & 
Bekki 2014; Zhukovska 2014; Bekki 2015a) as well as metal- 
licity (Inoue 2003; Asano et al. 2013a) in an attempt to bet¬ 
ter match the observed dust-metallicity relation. Following 
these prescriptions, for each gas cell we compute the local 
dust growth timescale 




(5) 


where p and T are the density and temperature of the gas 
cell, p ref and T ref are reference values for density and tem¬ 
perature in molecular clouds, and Tg ef is an overall normal¬ 
isation influenced by factors like atom-grain collision stick¬ 
ing efficiency and grain cross section (see Section 2.2 of Hi¬ 
rashita (2000) for a detailed derivation). We take p ref to 
be 1 H atom cm -3 and T ref = 20 K. This growth timescale 
is shortest in dense gas where collisions are more frequent 
than in the diffuse halo. Future efforts will need to more ac¬ 
curately model atom-grain collisions, taking into account a 
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possibly temperature-dependent collision sticking efficiency 
and a realistic grain size distribution (Weingartner V Draine 
2001; Li V Draine 2001). 

2.2.3 Grain Destruction 

Dust grains that have formed in the ISM can subse¬ 
quently be destroyed through a number of processes, in¬ 
cluding shocks from SN remnants (Seab V Shull 1983; Seab 
1987; Jones et al. 1994), thermal and nonthermal sputter¬ 
ing (Draine V Salpeter 1979b; Tielens et al. 1994; Caselli, 
Hartquist & Havnes 1997), and grain-grain collisions (Draine 
& Salpeter 1979a; Jones, Tielens & Hollenbach 1996). In par¬ 
ticular, shocks reduce the net efficiency of dust formation in 
SNe, expected to be the primary producers of dust at high 
redshift. The mass of grains destroyed in such a manner is 
thought to be proportional to the energy of the shocks (Mc¬ 
Kee 1989). 

Paralleling Equation (4), for every active cell we can 
estimate the local dust destruction rate for species i as 

d Mi 

,dust \ Mi ,dust /n\ 

d * /d T d 

where ra is a characteristic destruction timescale (McKee 
1989; Draine 1990; Jones & Nuth 2011) and again the dust 
mass is computed by summing its components originating 
from AGB stars, SNe la, and SNe II. While we could tie 
Td to the galaxy-wide SN rate (Hirashita & Ferrara 2002) 
or the radially-dependent gas surface density and SN rate 
(Calura, Pipino & Matteucci 2008), we instead estimate a 
dust destruction timescale on a cell-by-cell basis, given by 
Mg / 7 \ 

d e 7 M s ( 100 ) ’ 1 J 

where M g is the gas mass within a cell, e denotes the effi¬ 
ciency with which grains are destroyed in SN shocks, 7 is 
the local Type II SN rate, and M s (100) is the mass of gas 
shocked to at least 100 km s -1 (Dwek & Scalo 1980; Seab & 
Shull 1983; McKee 1989). We take e = 0.3, though there is 
a range of physically plausible values used in previous work. 
Applying the Sedov-Taylor solution to a homogeneous envi¬ 
ronment yields 

M s ( 100 ) = 6800.Esnii.si ( lnn p - 1 ) M ®- ( 8 ) 

\100 km s J 

where Esnii, 5 i is the energy released by a Type II SN in 
units of 10 51 erg, and v s ~ 100 km s -1 is the shock veloc¬ 
ity (McKee 1989). Equation 8 implicitly assumes SN shock 
expansion into a homogeneous medium of n — 0.13 cm -3 , 
corresponding to the star formation density threshold in 
our model. We employ a fixed ISM density when calculat¬ 
ing M s (100) because the detailed ISM multi-phase structure 
and, in particular, the ambient gas density around each SN 
are not resolved in our simulations. In any case, we caution 
that this destruction prescription neglects thermal sputter¬ 
ing and grain-grain collisions in the CGM and thus may 
overdeplete gas-phase metals in the diffuse halo. 

During every time step, the net dust growth rate in ev¬ 
ery active gas cell is computed by combining Equations (4) 
and ( 6 ), and this rate is used to update the local dust mass. 
When performing this update, we keep the relative propor¬ 
tions of dust mass originating from AGB stars, SNe la, and 
SNe II constant. 


Additionally, the amount of dust in the ISM is reduced 
when star particles are created. We assume that dust and 
gas-phase metals are distributed uniformly within each gas 
cell, so that each star particle’s metals come from gas-phase 
and dust sources in the same proportion as in the ISM. For 
example, if a star particle is created from a cell in which 25 
per cent of oxygen exists as dust and 75 per cent exists in the 
gas phase, then the amount of gas-phase oxygen lost by the 
cell is three times greater than the amount of oxygen dust 
lost. Other schemes to consume dust during star formation 
may be possible but more challenging to implement. 

2.2.4 Dust Effects on Cooling 

The galaxy formation model highlighted in Section 2.1 cal¬ 
culates gas cooling rates using contributions from primor¬ 
dial species, metals, and Compton cooling. In particular, 
the metal-line cooling rate is assumed to scale linearly with 
the gas-phase metallicity (see Equation (12) in Vogelsberger 
et al. 2013). The depletion of metals onto dust grains will 
decrease the gas-phase metallicity of the ISM and therefore 
reduce metal-line cooling. 

Observations suggest that dust cooling channels can ex¬ 
plain the formation of very metal poor stars (Caffau et al. 
2011; Klessen, Glover & Clark 2012), and modelling finds 
that cooling-induced fragmentation of dust impacts low- 
metallicity star formation (Schneider et al. 2006; Tsuribe 
& Omukai 2006; Dopcke et al. 2013). In numerical work, the 
temperature of dust grains can be computed via thermal 
equilibrium calculations accounting for atomic line emission, 
grain photoelectric effect heating, heating of dust through 
external radiation fields, and dust cooling through thermal 
emission (Krumholz, Leroy V McKee 2011). Dust may affect 
cooling rates computed in cosmological simulations using 
local ionising radiation (Cantalupo 2010; Gnedin V Hollon 
2012; Kannan et al. 2014). In hot gas with temperatures 
above 10 6 K, cooling from ion-grain collisions is expected 
to be strong (Ostriker V Silk 1973). For simplicity of our 
model, we neglect dust cooling channels in this first study. 

While our model does not currently implement any spe¬ 
cific dust cooling processes, metal-line cooling will be re¬ 
duced in comparison with previous simulations involving 
AREPO’s galaxy formation physics. Decreased cooling will 
in turn have a dynamical effect on galaxy formation, lower¬ 
ing the star formation rate (SFR). Future inclusion of dust 
cooling channels may weaken some of the dynamical effects 
seen in this work. 

2.2.5 Winds from Stellar Feedback 

We employ the non-local stellar feedback implementation 
from Vogelsberger et al. (2013) with some modification for 
dust. In this feedback prescription, gas cells in star-forming 
regions of the ISM are stochastically converted into wind 
particles, given a velocity 

v w -— ( 9 ) 

where r w is a dimensionless scale and ct^m is the local one¬ 
dimensional dark matter velocity dispersion, and allowed to 
move without hydrodynamic constraints. 

When wind particles recouple to the gas, they deposit 
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their metals, both gas-phase and dust, in the same propor¬ 
tion as the ISM from which they were launched. Thus, these 
stellar winds help carry metals and drive outflows away from 
dense regions of the ISM. In Section 4.2, we investigate the 
importance of these winds and their strength in distributing 
dust throughout a galaxy. 

3 COSMOLOGICAL RUNS 
3.1 Initial Conditions 

The initial conditions for our runs consist of the Aquarius 
suite of haloes, previously used for high-resolution cosmo¬ 
logical studies of Milky Way-sized structure (Springel et al. 
2008). All runs adopt the ACDM cosmological parameters 
T — 0.25, £2b — 0.04, Qa — 0.75, erg — 0.9, 
n s — 1, and Hubble constant Ho — 100/ikms -1 Mpc -1 = 
73 kms -1 Mpc -1 , which were employed in the Millennium 
and Millennium-II suites of simulations (Springel et al. 2005; 
Boylan-Kolchin et al. 2009). Though recent results from 
the Wilkinson Anisotropy Microwave Probe (Hinshaw et al. 
2013) and Planck satellite (Planck Collaboration et al. 2014) 
suggest, for example, a value of h < 0.73, we choose to adopt 
these parameters in part to remain consistent with recent hy- 
drodynamical simulations also utilising the Aquarius suite 
of haloes in AREPO (Marinacci, Pakmor & Springel 2014). 

The Aquarius haloes are labelled with the letters A 
through H, and each is simulated in a periodic box of side 
length 100 h -1 Mpc. In Aquarius nomenclature, our main 
simulations were performed using level 5 initial conditions. 
In addition, we simulated the Aquarius C halo using higher 
and lower resolution initial conditions, levels 4 and 6, re¬ 
spectively, in order to study the convergence properties of 
the fiducial dust model. We have found our dust model to 
be well-converged and use the highest resolution initial con¬ 
ditions for some visualisations in this work. The Aquarius C 
halo, which we also use to study variations in dust and feed¬ 
back physics, was adopted in the Aquila comparison project 
(Scannapieco et al. 2012) to analyse the results of a wide 
range of cosmological hydrodynamical codes. This halo has 
a fairly quiescent merger history, especially at low redshift 
(Wang et al. 2011). We also note that Marinacci, Pakmor V 
Springel (2014) applied this galaxy formation model without 
dust to the Aquarius haloes and found robust convergence 
over these same resolution levels. 

The same gravitational softening length was used for 
gas, dark matter, and stars, and it was kept constant in co¬ 
moving units until z = 1, with a maximum value in physical 
units of 680 pc for the high-resolution region using level 5 
initial conditions. The maximum values for level 4 and 6 
initial conditions were factors of two lower and higher, re¬ 
spectively. The z — 1 softening length was then used down 
to z — 0. We employ the SUBFIND algorithm (Springel et al. 
2001) for determining gravitationally-bound structure and 
substructure. 

In Table 1, we provide basic statistics on these haloes, 
computed at z — 0 using the fiducial model described below. 
We provide two computations of the gas and dust masses: 
M gas and Mdust, which are summed over all cells in the 
halo, and M gaS) disc and Mdust,disc, which include contribu¬ 
tions only from dense ISM gas cells and estimate the gas 


and dust content of the galactic disc. We isolate ISM gas 
by filtering cells according to temperature T and density p 
using the relation 

■° S (H) < 6 + °- 25 l0g ( 10.° [mA t pc- 3 ] ) ’ (10) 

which has been shown to remove cells in the hot halo (see 
Equation (1) in Torrey et al. 2012). Additionally, Table 1 
indicates the number and mass resolution of gas and dark 
matter cells in each halo. 

3.2 Fiducial Parameters 

In Table 2, we present the set of parameter values that com¬ 
prise our fiducial dust model. The fiducial feedback param¬ 
eters are similar to those used in Vogelsberger et al. (2013), 
with Esnii, 5 i = 1-09 and r w = 3.0, and we did not retune 
the feedback model after including dust. Given that dust 
depletion will affect cooling rates and star formation, future 
feedback modifications may be required. This highlights the 
need to include dust in detailed galaxy formation models. 

The fiducial dust condensation efficiencies follow those 
from Dwek (1998) and assume that a certain fraction of 
ejecta from an AGB star or SN exists as dust, with the re¬ 
mainder occupying the gas phase. These efficiencies are al¬ 
lowed to vary among species, since AGB stars and SNe are 
thought to produce dust of differing compositions (Ferrarotti 
& Gail 2006; Zhukovska, Gail V Trieloff 2008). For exam¬ 
ple, the high and low values of £^ GB ’ C// ° > 1 and respec¬ 
tively, are motivated by the idea that carbon-rich AGB stars 
are the dominant producers of carbonaceous grains. We cau¬ 
tion that the condensation efficiencies used for SNe la are 
high given recent observations and modelling work that sug¬ 
gest dust forms inefficiently in SNe la (Nozawa et al. 2011; 
Gomez et al. 2012). However, as shown in Section 4, SNe 
la produce only about 10 per cent of the dust in a galaxy, 
and so our results are not very sensitive to the condensation 
efficiencies for SNe la. 

While condensation efficiencies of the form described in 
Section 2.2.1 have also been adopted in more recent simula¬ 
tions (Calura, Pipino & Matteucci 2008; Bekki 2013), stel¬ 
lar models have begun to analyse the composition of ejected 
dust as a function of a star’s mass and metallicity (Ferrarotti 
& Gail 2006; Bianchi V Schneider 2007; Zhukovska, Gail 
V Trieloff 2008; Nanni et al. 2013; Schneider et al. 2014). 
For example, low-metallicity AGB stars are thought to pro¬ 
duce very low amounts of silicate grains (Ferrarotti V Gail 
2006). Also, some of these recent models track dust grains 
for SiC, AI 2 O 3 , or other compounds that do not map neatly 
onto our chemical evolution model. While we do not adopt 
condensation efficiencies that vary with mass or metallic¬ 
ity, we can still compare the dust masses predicted by our 
stellar yield tables and condensation efficiencies with those 
from more detailed stellar models. In Table 3, we calculate 
the total dust masses condensed in our model for various 
choices of stellar mass and metallicity. We note that AGB 
stars produce total dust masses around 10 -2 —10 -3 M 0 with 
no strong metallicity dependence, fairly consistent with stel¬ 
lar models (Zhukovska, Gail V Trieloff 2008; Nanni et al. 
2013). 

We can also compare the results of Table 3 to pre¬ 
dictions of dust masses condensed in SNe II, although we 
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Run 

■R 200 

[kpc] 

M gas 

[10 10 Mg] 

Ml 

[10 10 M 0 ] 

^Vlust 
[10 s M 0 ] 

■^gas,disc 

[10 10 M 0 ] 

^Viust,disc 

[10 s M 0 ] 

Vg as 

Ahm 

Ulgas 

[1O 5 M 0 ] 

771dm 
[10 5 M 0 ] 

A5:FI 

237.5 

11.62 

182.29 

5.61 

4.52 

3.79 

211527 

690478 

5.03 

26.40 

B5:FI 

183.0 

4.22 

78.18 

3.48 

1.37 

1.82 

119386 

518981 

3.35 

17.59 

C5:FI 

233.5 

11.11 

175.93 

5.98 

3.47 

3.47 

253836 

814834 

4.11 

21.59 

D5:FI 

240.7 

15.18 

195.56 

12.20 

6.12 

7.66 

319074 

846419 

4.40 

23.10 

E5:FI 

206.7 

6.14 

114.16 

5.35 

1.09 

1.00 

179039 

652270 

3.33 

17.50 

F5:FI 

208.0 

9.26 

113.67 

9.48 

3.58 

5.55 

375068 

942365 

2.30 

12.06 

G5:FI 

201.9 

9.94 

95.99 

7.78 

5.20 

5.74 

317798 

769854 

2.83 

14.88 

H5:FI 

180.1 

2.61 

76.24 

1.67 

0.20 

0.34 

87369 

588050 

2.96 

15.56 

C4:FI 

232.5 

6.82 

159.19 

3.89 

1.44 

1.76 

1265814 

5898234 

0.51 

2.70 

C6:FI 

235.6 

13.92 

179.84 

6.91 

4.97 

4.54 

39592 

104118 

32.90 

172.73 


Table 1. Basic data on all eight Aquarius haloes simulated using level 5 initial conditions and the fiducial dust model detailed in 
Section 3.2. Additionally, we show the same statistics for the Aquarius C halo simulated at levels 4 and 6. To keep consistent with the 
naming convention adopted later in this work, runs are referred to by their halo letter (A through H), resolution level (4, 5, or 6), and 
the suffix “FI” to denote that the fiducial dust physics were used. Here, the virial radius R 200 is the radius about the halo’s potential 
minimum enclosing a density 200 times greater than the critical density. We list the gas, dark matter, and dust masses using all cells 
within the halo, denoted M gas , M & m , and Md us t, respectively. Here, M gas ^- lsc and Md us t,disc are the estimated gas and dust masses in 
the galactic disc, computed using only dense ISM gas cells as determined by the temperature-density phase space cut in Equation (10). 
Finally, we indicate the number of gas and dark matter cells in the halo and their respective mass resolutions. 


Variable Fiducial Value 


Description 


.AGB,C/O > 1 

0.0 

for i 


1.0 

for i 

^AGB,C/0 < 1 

0.0 

for i 


0.8 

for i 

<5f N 

0.0 

for i 


0.5 

for i 


0.8 

for i = 

T-ref 

g 

0.2 


Esnii,5i 

1.09 


H, He, N, O, Ne, Mg, Si, Fe, 
C 

H, He, N, C, Ne, 

O, Mg, Si, Fe 
H, He, N, Ne, 

C, 

O, Mg, Si, Fe 


dust condensation efficiency for species i 

in AGB stars with C/O > 1 in ejecta 

dust condensation efficiency for species i 

in AGB stars with C/O < 1 in ejecta 

dust condensation efficiency for species i in SNe 


reference dust growth timescale, in units of [Gyr] 
available energy per SN II, in units of [10 51 erg] 


Table 2. Fiducial dust model parameters. The dust condensation efficiencies regulate the fraction of metal ejecta for a chemical species 
that condenses into dust, as outlined in Section 2.2.1. The reference growth timescale r g ei influences the rate at which gas-phase metals 
deplete onto dust grains, and Esnii, 51 controls the energy per SN II, which in turn affects grain destruction in SN shocks. 


caution that this subject is debated in the literature. For 
SNe II with masses of 12 M© and 20 M©, we calculate con¬ 
densed dust masses of roughly 0.4 M© and 0.8 M©, respec¬ 
tively. In comparison, Todini <V Ferrara (2001) predicted 
about 0.3 M© and 0.5 M© of dust for these choices of stel¬ 
lar mass. Our results are consistent with the estimate that 
between 2 per cent and 5 per cent of a progenitor’s mass 
becomes dust in SNe II in the range 13 M© ^ M ^ 40 M© 
(Nozawa et al. 2003). However, there is some tension with 
observations: SN 2003gd, for example, was determined to 
have no more than 0.02 M© of dust (Sugerman et al. 2006; 
Meikle et al. 2007), and other SNe have failed to produce 
the amount of dust expected from models (Rho et al. 2008, 
2009). As stellar models begin to implement more complex 
chemical reaction networks (e.g. Cherchneff 2010) and more 
observations are gathered, better constraints can be placed 
on the dust masses condensed in core-collapse SNe. 

In the ISM, refractory grains are estimated to have de¬ 
struction timescales of 10 8 — 10 9 yr, lower than the typi¬ 
cal injection time of dust from stellar sources (Barlow 1978; 
Draine Salpeter 1979a; Dwek Scalo 1980; McKee 1989). 
Mechanisms that allow metals to recondense in the ISM 
and counter the destruction-injection imbalance are needed 


to regulate gas-phase abundances. The growth timescale r g 
characterises dust growth in molecular clouds and the asso¬ 
ciated probability that a gas-phase metal atom will collide 
and stick to an existing grain (Draine 1990; Hirashita 2000). 
Just as we tie the dust destruction timescale in each cell to 
its local SN rate, we adopt a growth timescale dependent on 
local gas density and temperature that enables dust to grow 
more quickly in dense ISM gas. There is evidence that dust 
growth is particularly dominant over stellar dust production 
in galaxies above a certain critical metallicity (Inoue 2011; 
Asano et al. 2013a; Mancini et al. 2015), and other mod¬ 
els have begun to explore the effect of density, temperature, 
and metallicity variations on local growth timescales (In¬ 
oue 2003; Bekki 2015a). While models that employ variable 
growth timescales sometimes assume a characteristic grain 
size, grain density, and atom collision sticking efficiency, in¬ 
creasing the number of estimated parameters, in this work 
we fixed such parameters at typical Galactic values in a man¬ 
ner similar to Equation (12) in Hirashita (2000) in order to 
capture the essential density and temperature dependence 
in the growth timescale. The resulting estimated reference 
timescale r g ef is shown in Table 2. 

Table 4 lists the range of model variations that we con- 


© 2016 RAS, MNRAS 000 , 1-28 









8 R. McKinnon, P. Torrey, and M. Vogelsberger 


M [Mg] 

M dust (Z = 0.004) [M 0 ] 

M dust (Z = 0.008) [M e ] 

M dust (Z = 0.02) [M 0 ] 

2 

0.065 

0.025 

0.009 

4 

0.003 

0.006 

0.015 

6 

0.006 

0.011 

0.029 

12 

0.392 

0.381 

0.365 

20 

0.780 

0.781 

0.780 


Table 3. Sample of condensed dust masses (Mdust) as a function of initial stellar mass (M) and metallicity (Z) using the stellar yields 
and dust condensation formulae in Section 2. We assume stars have scaled solar abundances. In our simulations, the transition between 
AGB stars and SNe II occurs at M = 8Mq. The yields for AGB stars are presented in Karakas (2010) at Z = 0.019 but listed here as 
Z = 0.02 for easy comparison with SNe yields at Z — 0.02. 


Name 

Abbreviation 

Physics 

fiducial 

FI 

fiducial dust parameters from Table 2, k w = 3.0, includes AGN feedback 

no feedback 

NF 

no stellar or AGN feedback 

AGN feedback 

AF 

only AGN feedback 

slow winds 

SW 

/■c-w 1.85 

fast winds 

FW 

G 

II 

£ 

se 

constant growth timescale 

CG 

uses growth timescale of r g = 0.2 Gyr, ignores local density and temperature 

constant destruction timescale 

CD 

uses destruction timescale of r d = 0.5 Gyr, ignores local SN rate 

low production 

LP 

all condensation efficiencies from Table 2 halved 

no growth 

NG 

no growth mechanism 

no destruction 

ND 

no destruction mechanism 

production only 

PO 

no growth or destruction mechanisms 

without dust 

WD 

no dust tracking, all metals exist in gas phase 


Table 4. Description of different runs with varying feedback and dust model parameters. The first and second columns characterise the 
nature of the model variation, and the third column specifies the exact physics changes. In particular, the no feedback (NF) and AGN 
feedback (AF) models ignore all stellar feedback and winds. The slow winds (SW) and fast winds (FW) runs vary the dimensionless 
parameter k, w responsible for scaling the velocity of wind particles. The constant destruction timescale (CD) variation adopts a constant 
value of r d used in Equation (6), mimicking previous one-zone models that ignore local SNe when destroying grains. The low production 
(LP) run adopts condensation efficiencies that are half of those of the fiducial model in an attempt to assess the importance of the 
grain formation process. The without dust (WD) model uses fiducial feedback parameters but does not track dust. It excludes all of the 
changes in Section 2.2. 


sider in addition to the fiducial setting. These variations are 
divided into two categories: those that apply the fiducial 
dust model and vary the underlying feedback mechanisms, 
and those that adopt the fiducial feedback settings and ex¬ 
plore the impact of the various dust processes. We also ex¬ 
plore the fiducial feedback model without dust tracking. In 
this last model, all metals are assumed to exist in the gas 
phase and contribute to metal-line cooling. This model with¬ 
out dust is expected to produce larger star formation rates 
and gas-phase metallicities. 

These feedback and dust model variations, together 
with the range of Aquarius haloes used as initial conditions, 
provide a suitable starting point for understanding the im¬ 
pact of dust on cosmological galaxy formation simulations. 

In the remainder of this paper, we refer to our runs us¬ 
ing four-character identifiers: the first refers to the Aquarius 
halo chosen (i.e. one of A through H), the second denotes the 
resolution level of the initial conditions (i.e. 4, 5, or 6), and 
the last two indicate the underlying physical model using 
the abbreviations in Table 4. For example, C5:ND is short¬ 
hand for the simulation of the Aquarius C halo using level 
5 initial conditions with no dust destruction mechanism. 


4 RESULTS 

4.1 Distribution of Dust in the Fiducial Model 

We first use our highest resolution simulation of the Aquar¬ 
ius C halo with fiducial dust and feedback physics to anal¬ 
yse the distribution of gas-phase metals and dust within 
the central halo and surrounding CGM. This is motivated 
by observations of significant amounts of dust in the CGM 
(Menard et al. 2010; Peeples et al. 2014; Peek, Menard & 
Corrales 2015) and the fact that gas-phase metals may con¬ 
dense into dust in regions of low star formation. In Figure 1, 
we show surface densities of baryons (E), gas-phase metals 
(Eg as _phase metal), and dust (E dus t) for face-on and edge-on 
views of the galactic disc at z = 2, 1, and 0. We also dis¬ 
play the projected dust-to-metal ratio (Z dust /Z) at these 
redshifts, where Z du st denotes the mass fraction of dust and 
Z is the usual metallicity, including gas-phase metals and 
dust. This halo does not undergo any major mergers below 
z — 6 (Wang et al. 2011), and so these images capture the 
quiescent formation of a disc of diameter roughly 15kpc. 

At z = 2, dust is most concentrated in the galactic 
center, with a surface density roughly an order of magnitude 
larger than that outside the disc. The surface densities of 
gas-phase metals and dust evolve in a similar fashion from 
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XI [Mq pc ] ^gas-phase me tal [^0 PC ] list. [Mq P^ ] ^dust/-^ 

Figure 1 . Surface densities of baryons, gas-phase metals, and dust (first, second, and third columns, respectively) as well as corresponding 
dust-to-metal ratio (fourth column) for the Aquarius C halo at z = 2, 1, and 0 using level 4 initial conditions with fiducial feedback and 
dust physics. The two images in each column at fixed redshift represent face-on (top) and edge-on (bottom) projections. All distances are 
given in physical units. The scale bar in the first row indicates 25kpc. Projections were performed in a cube of side length 150 1 kpc 

centered on the potential minimum. Movies of this simulation are available through http://www.mit.edU/~ryanmck/#research. 
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Figure 2. Zoomed-in projections of the dust-to-metal ratio using the face-on (top) and edge-on (bottom) views of the Aquarius C halo 
at z = 2, 1, and 0 presented in Figure 1. These plots capture the inner disc region, with the scale bar in the upper right indicating lOkpc 
in physical units. A cube of side length 50 1 kpc was used for the projection volume. 


z — 2 to z — 1 and largely trace the distribution of baryons. 
By z — 1, Edust > 10 -2 M© pc -2 extends out to around 
25 kpc as dust grows on the edges of the galactic disc. The 
surface density of gas-phase metals in the CGM is largely 
unchanged from z — 2 to z — 1, with a dust-to-metal ratio 
below 0.1. The dust-to-metal ratio is largest near the galactic 
center and several times higher than the typical value away 
from the disc. By z = 0, there is less variation between 
gas-phase metals and dust. The central disc, roughly 15 kpc 
across, has Edust > 10 -1 M©pc -2 , and the dust-to-metal 
ratio remains highest in the star-forming central core where 
dust is produced. We caution that the absence of thermal 
sputtering in the diffuse halo may artificially cause dust-to- 
metal ratios to rise in the CGM. 


These results indicate that dust depletion can affect 
the distribution of gas-phase metals, especially in regions 
where high star formation activity produces the most dust. 
To highlight the dust evolution in the inner galactic disc, in 
Figure 2 we display zoomed-in projections of the dust-to- 
metal ratio for the face-on images from Figure 1. At each of 
z = 2, 1, and 0, the dust-to-metal ratio is largest near the 
galactic center. However, by z — 0, the spatial variation in 
dust-to-metal ratio that is prominent at z — 1 diminishes, 
although there are small pockets of low dust-to-metal ratio 
in the outer disc. In conjunction with Figure 1, this suggests 
that observations of galaxies should see dust-to-metal ratios 
that are highest near galactic cores and vary with redshift. 


4.2 Impact of Feedback 

While some previous cosmological simulations treating dust 
have investigated the effect of galactic winds (Zu et al. 2011) 
and SN feedback (Bekki 2015a), there has been no compre¬ 
hensive analysis of feedback physics on dust evolution. In 
Figure 3, we show the total dust, metal, gas, and stellar 
masses as a function of redshift for each of the feedback 
variations detailed in Table 4 applied to the Aquarius C 
galaxy. We also plot the SFR and gas-phase metallicity, two 
dynamical quantities that will be impacted by the presence 
of dust. 

Until z ~ 4, the no feedback and AGN-only feedback 
models are very similar, producing stellar masses several 
times larger than those seen in the fiducial, slow winds, 
and fast winds models with stellar feedback enabled. This is 
consistent with the SFR plot, which shows that the no feed¬ 
back and AGN-only feedback runs have similar behaviour at 
high redshift. Ah runs yielded dust masses at high redshift 
more than an order of magnitude smaller than the several 
10 7 M© and above of dust observed for SDSS J1148+5251 
and A1689-zDl at z — 6.4 and higher (Valiante et al. 2009, 
2011; Watson et al. 2015). This is not too surprising since 
simulations of haloes larger than those in the Aquarius suite 
seem to be required to produce such dusty systems at high 
redshift. 

At low redshift, AGN feedback strongly suppresses the 
SFR in the Aquarius C galaxy. The reduction in SFR pro¬ 
vides more opportunity for gas-phase metals in the ISM to 
condense into dust, leading to a decline in gas-phase metal- 
licity relative to the no feedback model. The runs with stellar 
feedback see a similar suppression of gas-phase metallicity 
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log(l + z) 


log(l + z) 


Figure 3. Clockwise from top left, a comparison of dust mass, metal mass, stellar mass, gas-phase metallicity, star formation rate, and 
gas mass as a function of redshift for the Aquarius C galaxy using various feedback parameters. These quantities include contributions 
only from the dense galactic disc and not the hot halo. In each plot, the top panel shows results on an absolute scale, while the ratio of 
quantities relative to the fiducial run (denoted /pi) is shown in the bottom panel. This convention is also adopted in subsequent figures. 
The metal mass includes both gas-phase and dust contributions. The fast winds (FW) and no feedback (NF) runs produce dust masses 
at 2 : = 0 roughly half that produced in the fiducial (FI) model. The fast winds run sees the strongest decline in SFR and gas-phase 
metallicity towards z = 0. 
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Figure 4. Comparison of dust-to-gas ratio (T>; left) and dust-to-metal ratio ( Z^ us i/Z ; right) for the range of feedback models run on 
the Aquarius C halo. Quantities correspond to the dense galactic disc, as in Figure 3. The fiducial (FI), slow winds (SW), and fast winds 
(FW) runs with stellar feedback yielded similar dust-to-gas ratios at all redshifts. There is significant scatter in the dust-to-metal ratio, 
though the feedback run with the largest star formation rate and gas-phase metallicity at 2 : = 0 displayed the smallest dust-to-metal 
ratio. This is consistent with our SN-driven dust destruction mechanism. 


below z ~ 1. This effect is most pronounced for the fast 
winds model, which is most efficient at driving gas-phase 
metals away from star-forming regions and has the lowest 
stellar mass. In the absence of SN activity, these gas-phase 
metals can accrete onto dust grains more rapidly. The fast 
winds run sees a drop in gas-phase metallicity of roughly 
0.7 dex from z — 1 to z = 0. We note that the gas-phase 
metallicities shown in Figure 3 correspond to the dense 
galactic disc and that the hot halo sees less metal enrich¬ 
ment. The fiducial and slow winds models do not experience 
such severe declines in gas-phase metallicity and are more 
similar to the AGN-only feedback run at z — 0 than the 
extreme no feedback model. 

Figure 4 shows the evolution of the dust-to-gas ratio D 
and dust-to-metal ratio Z^ust/Z for the same feedback vari¬ 
ations. The fiducial model for this Milky Way-sized galaxy 
yields D « 10 -2 , consistent with estimates of the Galac¬ 
tic value (Gilmore, Wyse & Kuijken 1989; Sodroski et al. 
1997; Zubko, Dwek & Arendt 2004). The feedback variations 
also yield dust-to-gas ratios within a factor of several of the 
fiducial value. The dust-to-gas ratio increases by roughly 
1 dex from z = 2 to z = 0 for the fiducial, slow winds, and 
fast winds runs. Above, we noted that the fast winds model 
promotes gas-phase depletion of metals, and so one might 
expect a high dust-to-gas ratio. However, because star for¬ 
mation is so strongly suppressed in this model, the total 
amount of dust is roughly half that of the fiducial run. As 
a result, the fast winds model has a dust-to-gas ratio com¬ 
parable to the fiducial value. On the other hand, the fast 
winds z = 0 dust-to-metal ratio is the largest of all models, 
at nearly 0.75. This is slightly below the fiducial model’s 


z = 0 dust-to-metal ratio of near 0.65, which is still above 
but closer to estimates of the Galactic value. The no feed¬ 
back run with more star formation was more effective at 
regulating the dust-to-metal ratio and returning dust to the 
gas phase. 

The dust-to-metal ratio is expected to increase signifi¬ 
cantly at low redshift: previous modelling predicts Z^ust/Z 
at z = 0.5, 1, and 2 to be 50, 30, and 20 per cent of the 
z — 0 value (Inoue 2003). While the feedback variations do 
not reproduce these precise numbers, there is a noticeable 
increase towards low redshift. This increase is present across 
ah feedback variations, even though the z — 0 dust-to-metal 
ratios range from 0.2 for the no feedback model to over 0.7 
for the fast winds model. The presence of feedback affects 
the normalisation of the dust-to-metal ratio but does not 
strongly alter the shape of its evolution. 

While the dust-to-metal ratio was strongly affected by 
variations in feedback, other quantities, like surface density 
of dust projected onto the galactic plane, are less sensitive. 
In Figure 5, we show the dust surface density Hdust versus 
radial distance at z = 0 for each feedback model. Outside of 
the galactic disc, we compare with the Edust oc r -0 ' 8 scaling 
observed in SDSS data (Menard et al. 2010). Additionally, 
since the Aquarius suite of haloes forms Milky Way-sized 
galaxies, we overlay recent observations of the dust surface 
density in the inner disc of M31, which has an estimated 
dust mass of 5.4 x 10 7 M© (Draine et al. 2014). In con¬ 
trast to previous figures, ah feedback variations yield simi¬ 
lar results, with the surface density profiles peaking around 
Sdust ~ 10 -1 M©pc -2 for z — 0. These surface densities 
display moderate increases from z = 1 to z = 0, consistent 
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Figure 5. Surface density of dust (Edust) as a function of radial 
distance from the Aquarius C halo’s spin axis at z = 0 for our 
set of feedback variations. The radial range extends to R 200 , as 
defined in Table 1, and the bottom subpanel shows surface den¬ 
sity relative to the fiducial result. The gray dashed line indicates 
the Edust c* r -0-8 scaling observed in SDSS data over cosmologi¬ 
cal distances (Menard et al. 2010). The amplitude of this scaling 
has been chosen to align with simulated data from 25 kpc out 
to roughly 100 kpc. Squares indicate measurements of M31’s pro¬ 
jected dust profile (Draine et al. 2014). Triangles denote the M31 
data scaled by a factor of two to model a Milky Way dust mass 
that may be higher than that of M31. 


with the growth in dust mass from Figure 3. All feedback 
runs experience a drop in dust surface density of roughly 
one order of magnitude or more when leaving the galac¬ 
tic disc, with normalisations and shapes in broad agree¬ 
ment with M31 data. These models also all predict Edust 
scalings out to more than 100 kpc that are fairly consistent 
with SDSS observations. However, cosmological simulations 
of larger galaxy populations will be needed to confirm this 
trend on the scales detailed in Menard et al. (2010). Such 
simulations will also be able to investigate predictions about 
the dust content of the intergalactic medium (Aguirre 1999; 
Inoue V Kamaya 2003; Petrie et al. 2006; Dijkstra V Loeb 
2009; Corrales & Paerels 2012). 

Figure 6 displays the dust-metallicity relation for these 
feedback runs, plotting dust-to-gas ratio versus gas-phase 
metallicity at z = 2, 1, and 0. We also display observations 
from Kennicutt et al. (2011) and Galametz et al. (2011), 
which cover roughly a 1.5 dex range in oxygen abundance. To 
facilitate comparison with these data, we compute the dust- 
to-gas ratio and gas-phase metallicity using only dense gas in 
the galactic disc, as determined by the temperature-density 
phase space cut in Equation (10). We use this temperature- 
density cut in all subsequent plots involving comparison to 
observational data. 

Variations in feedback affect both the gas-phase metal¬ 
licity and the dust-to-gas ratio, particularly at low redshift. 


Figure 6. Dust-to-gas ratio ( D ) versus gas-phase metallicity for 
varying feedback runs involving the Aquarius C halo, plotted at 
z = 2,1, and 0. Smaller circles denote higher redshift. The dust- 
to-gas ratio and gas-phase metallicity have been computed using 
only dense gas cells in the galactic disc. For comparison, obser¬ 
vational data from Kennicutt et al. (2011) and Galametz et al. 
(2011) (as compiled in Remy-Ruyer et al. 2014) are provided. All 
runs except that with no feedback experience a drop in gas-phase 
metallicity from z = 1 to z = 0. 


All models except the no feedback model see a decline in 
gas-phase metallicity from z = 1 to z = 0, with the effect 
most pronounced for the fast winds run. The fiducial model 
is most similar to the runs with fast and slow winds, con¬ 
sistent with the results in previous figures. Although these 
stellar feedback models lie above the dust-metallicity rela¬ 
tion at z = 0, the fits at z — 1 and z — 2 are much more 
reasonable. The no feedback model has a high gas-phase 
metallicity at all redshifts and thus lies slightly below the 
dust-metallicity observations. These results indicate that the 
inclusion of feedback in our dust model can affect the dust- 
metallicity evolution of galaxies and help limit high-redshift 
metallicities. 

Because heavy elements trapped in dust grains do not 
contribute to metal-line cooling in our model, we expect star 
formation to decrease as gas-phase metals accrete onto dust. 
This qualitative effect is shown in Figure 7, which compares 
the projected baryon and gas-phase metal densities for our 
fiducial model as well as a model with all dust physics dis¬ 
abled. In this latter model, all metals outside of stars are 
considered to be in the gas phase. The run without dust pro¬ 
duces a slightly larger galactic disc at z = 0, while the run 
with dust slightly reduces the density of gas-phase metals in 
the CGM. While we caution that our fiducial dust model ig¬ 
nores dust sputtering and may overdeplete gas-phase metals 
in the halo, these results suggest that the presence of dust 
affects how galaxies evolve in simulations, even if feedback 
settings are unaltered. 

To analyse this effect quantitatively, in Figure 8 we plot 
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Figure 7. Face-on and edge-on projections of total surface den¬ 
sity (left column) and gas-phase metal surface density (right col¬ 
umn) at z = 0 for simulations of the Aquarius C halo with fidu¬ 
cial dust (FI; top half) and without dust (WD; bottom half). 
The scale bar in the first row marks 25kpc, and the projection 
volume is the same as that used in Figure 1. Not shown are met¬ 
als occupied by dust, which exist only in the fiducial simulation. 
Compared to the simulation without dust, the fiducial dust run 
results in a smaller disc and a decreased gas-phase metal surface 
density beyond 15 kpc from the galactic center. 
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Figure 8. A comparison of star formation rate (SFR; top) and 
gas-phase metallicity (bottom) computed for the Aquarius C 
galaxy as a function of redshift in the fiducial (FI) and without 
dust (WD) models. The FI run has decreased metal-line cooling, 
as gas-phase metals deplete onto dust grains. The SFRs largely 
agree, but at z = 0, the difference in gas-phase metallicities be¬ 
tween the FI and WD runs is roughly 0.5 dex. 


the SFR and gas-phase metallicity versus redshift for our 
fiducial model as well as our model without dust. Differences 
in these models are most noticeable beginning at z « 1. 
While both models produce similar SFRs, the deviation in 
gas-phase metallicity from z ~ 1 to z — 0 is more pro¬ 
nounced. Without dust, the gas-phase metallicity increases 
and plateaus, whereas the model with dust sees a decline of 
roughly 0.5 dex. At low redshift, even the fiducial feedback 
settings are sensitive to whether dust is included or not. 

In the analysis above, the fiducial feedback model pro¬ 
duced dust-to-gas and dust-to-metal ratios in rough agree¬ 
ment with the Galactic values. While its location on the 
observed dust-metallicity relation was in tension with ob¬ 
servations at z — 0, the fiducial run at higher redshift 


yielded more accurate results. The no feedback model, al¬ 
ready known to overproduce high stellar mass galaxies (Vo¬ 
gelsberger et al. 2013), has the highest SFR and lowest dust- 
to-metal ratio of the various models. Fast winds strongly 
suppress star formation and dust production, carrying met¬ 
als to galactic regions where dust destruction is less effective. 
This resulted in a poor z = 0 dust-metallicity fit, owing to 
low redshift gas-phase depletion. While the models with slow 
winds and AGN-only feedback do not suffer such extreme 
problems, they offer no advantage over the fiducial feedback 
model. Previous studies have shown that in the absence of 
winds, the dust-to-metal ratio required to match observed 
intergalactic reddening is unphysical (Zu et al. 2011). We 
note that variations in stellar wind feedback affect our dust- 
metallicity relation more strongly than in Bekki (2015a), 
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Figure 9. Dust mass as a function of redshift for a variety of 
dust models applied to the Aquarius C galaxy. This only includes 
contributions from the inner disc and not the hot hao. The no 
feedback (NF) run is shown for comparison. At z = 0, the dust 
model without a growth mechanism (NG) yields a final dust mass 
more than an order of magnitude below those obtained by models 
with growth enabled. The difference between a destruction mech¬ 
anism based on local SNe (FI) and a uniform timescale (CD) is 
minor at z = 0, but at high redshift the CD run is more simi¬ 
lar to the run without any destruction (ND). The model without 
any growth or destruction (PO) sees a dust mass at low redshift 
reduced by roughly a factor of two compared to the ND run with 
growth. The fiducial model was most similar to the low produc¬ 
tion (LP) model with smaller condensation efficiencies. 

though in both cases stronger feedback suppresses the dust- 
to-gas ratio and gas-phase metallicity. In the remainder of 
this paper, we adopt the fiducial feedback parameters. 

4.3 Variations of the Dust Model 

Motivated by observations of high-redshift dusty galaxies, 
previous one-zone dust models have studied the importance 
of dust growth in the ISM relative to contributions from 
stellar sources. Dust growth appears to be dominant in 
galaxies above a critical metallicity, and galaxies with a 
shorter star formation timescale reach this critical metal¬ 
licity more quickly (Inoue 2011; Asano et al. 2013a). The 
carbon-to-silicate dust mass ratio changes significantly as 
dust growth overtakes stellar injection, and dust growth in 
the ISM can be important even in dwarf galaxies (Zhukovska 
2014). While we do not investigate it in this paper, we ac¬ 
knowledge that variations in the IMF can affect dust model 
results, particularly in low mass galaxies (Gall, Andersen & 
Hjorth 2011a). Top-heavy IMFs may also strengthen dust 
production from SNe, despite leading to increased grain de¬ 
struction (Dunne et al. 2011). It is therefore natural to con¬ 
sider how the various components of our dust model impact 
galaxy evolution. This will help to determine which aspects 


of dust physics are most important and worthy of more de¬ 
tailed modelling in future work. 

In Figure 9, we plot the redshift evolution of dust mass 
in the Aquarius C galaxy for each of the dust model vari¬ 
ations listed in Table 4. We also show how each variation 
compares to the fiducial result. Because dust was shown to 
have a dynamical effect on star formation in Section 4.2, 
the total metal masses for each run show minor differences. 
The dust mass evolution of the fiducial model was most 
similar to that of the model with decreased dust conden¬ 
sation efficiencies, enabling more mass to be ejected in the 
gas phase. In fact, their z = 0 dust masses were almost 
identical, while at z ~ 2.5 the fiducial model produced 
about twice as much dust. This perhaps reduces the effect 
of large uncertainties in dust condensation efficiencies, es¬ 
pecially for results at low redshift. For modelling at high 
redshift, metallicity-dependent condensation efficiencies (as 
in the one-zone model of Zhukovska, Gail & Trieloff 2008) 
may have more of an impact. 

In contrast, the model without dust growth deviated 
significantly from the fiducial run, and the dust mass in¬ 
creased by only a factor of four from z 4 to z = 0. The 
fiducial model and other models with growth enabled saw 
an increase in dust mass closer to one order of magnitude. 
This corroborates findings from earlier one-zone models and 
highlights the importance of dust growth at low redshift. 
We caution that more detailed growth mechanisms follow¬ 
ing one-zone models that take into account variations in 
dust grain sizes (Hirashita & Kuo 2011; Asano et al. 2013b; 
Nozawa et al. 2015) and their impact on grain collision rates 
will offer a better assessment of the role dust growth plays 
in the ISM. Before z « 5, the difference between a model 
employing only stellar dust production and a model with 
both stellar production and ISM dust growth is minor, while 
both models yield more dust than in the fiducial run with 
destruction. However, by z = 0, the model lacking growth 
and destruction sees a final dust mass approximately two 
times smaller than that obtained by the model with growth 
but without destruction. Grain growth becomes increasingly 
important towards low redshift. This is consistent with the 
hypothesis that galaxies pass through a critical metallicity 
at which dust growth in the ISM dominates stellar produc¬ 
tion. 

The run adopting a constant destruction timescale of 
Td = 0.5 Gyr for dust destruction, rather than the fiducial 
timescale computed locally using SN activity, behaves more 
like the model lacking any destruction whatsoever and over¬ 
produces dust at high redshift relative to the fiducial model. 
By z = 0, the constant destruction timescale run agrees 
well with the fiducial run. Dust models that globally adopt 
Td = 0.5 Gyr (e.g. the fiducial runs in Bekki 2015a) may be 
overestimating a galaxy’s dust content by a factor of several 
at high redshift. 

Figure 10 shows the dust-to-gas and dust-to-met al ra¬ 
tios for these dust model runs, using contributions from cells 
within the disc of the Aquarius C galaxy. Nearly all dust 
variations yield a z = 0 dust-to-gas ratio within a factor of 
two of the fiducial value of D « 10 -2 . At high redshift, mod¬ 
els without dust destruction display an increased dust-to-gas 
ratio. A dust model without growth results in a late time 
dust-to-gas ratio suppressed by over 1 dex relative to the 
fiducial run. Without growth, the dust-to-gas ratio shows 


© 2016 RAS, MNRAS 000, 1-28 








16 


R. McKinnon, P. Torrey, and M. Vogelsberger 


t [Gyr] 



t [Gyr] 



log(l + z) 


log(l + z) 


Figure 10. Comparison of dust-to-gas ratio ( D ; left) and dust-to-metal ratio (^dust/^5 right) for the range of dust models in Figure 9 
run on the Aquarius C galaxy. These ratios are computed using total masses within the dense galactic disc. The no feedback (NF) model 
is also shown for comparison. In both plots, the model lacking a dust destruction mechanism (ND) yielded the highest dust-to-gas and 
dust-to-metal ratios at essentially all redshifts. Conversely, the model with only dust production enabled (PO) results in dust-to-gas and 
dust-to-metal ratios at z = 0 roughly half those obtained by the fiducial model. The run without dust growth but with destruction (NG) 
more strongly suppresses both ratios. The NG and PO runs display very flat dust-to-metal ratio evolutions. 


little evolution from its value at z ~ 5. The effect of limiting 
dust growth is much more severe than that of limiting grain 
destruction. Modelling work motivated by observations of 
SINGS galaxies and high redshift quasars has similarly pre¬ 
dicted that dust growth and non-stellar dust production can 
have a pronounced impact on a galaxy’s dust content, while 
dust destruction is less influential (Draine 2009; Michalowski 
et al. 2010; Valiante et al. 2011; Mattsson, Andersen & 
Munkhammar 2012; Zafar & Watson 2013; Rowlands et al. 
2014). 

Models without growth and without destruction pro¬ 
duced the lowest and highest dust-to-metal ratios for z < 2, 
respectively. While the dust-to-metal ratio predicted by only 
allowing stellar production of dust is roughly 0.3 and phys¬ 
ically acceptable, it is nearly constant and displays little 
evolution. In contrast, both growth and destruction mech¬ 
anisms are needed to capture the increase in dust-to-metal 
ratio expected at low redshift (Inoue 2003) , as in our fiducial 
model and the run with low dust condensation efficiencies. 
The extent to which the dust-to-metal ratio varies within a 
galaxy and among different galaxies is unclear (Mattsson, 
Andersen & Munkhammar 2012), so it is difficult to assess 
whether our fiducial dust model yields a typical dust-to- 
metal ratio. In any case, from Figure 10 we can conclude that 
varying dust condensation efficiencies makes little difference 
to the dust-to-metal ratio evolution, and while adopting a 
constant destruction timescale instead of using local SNe- 
driven destruction leads to a higher dust-to-metal ratio at 
high redshift, the effect is nowhere near as severe as limiting 
dust growth or destruction. 


In Figure 11 we plot the dust-to-gas ratio versus dust- 
to-metal ratio at z — 2, 1, and 0 for every gas cell in the 
halo using the fiducial run, a model without dust growth, 
and a model without dust destruction. Given the impor¬ 
tance of metallicity in regulating the transition between 
stellar production-dominated dust evolution and ISM-led 
dust growth, we also compute the mean metallicity for each 
region of dust-to-gas and dust-to-metal ratios. While the 
Aquarius C halo is fairly quiescent over this redshift range, 
we see qualitative differences between these three dust mod¬ 
els with slight evolution over time. 

At all redshifts, the fiducial model contains cells across 
a wide range of dust-to-metal ratio. Regions of larger dust- 
to-gas ratio tended to be more metal-rich, with the peak 
dust-to-gas ratio and metallicity increasing by roughly 1 dex 
and 0.3 dex, respectively, from z — 2 to z — 0. In the fiducial 
model, low dust-to-gas ratios at z — 0 were concentrated in 
cells with low dust-to-metal ratio. While one-zone dust mod¬ 
els are incapable of producing such phase space plots, previ¬ 
ous smoothed-particle hydrodynamical simulations have ex¬ 
plored the relation between dust-to-gas and dust-to-metal 
ratios for fiducial dust parameters (Bekki 2015a). While our 
results reproduce the overall increase in dust-to-gas ratio 
from z = 2 to z = 0 and roughly 1 dex scatter in typical 
dust-to-gas ratio at z = 0, previous work finds nearly all 
gas cells concentrated in a narrow band of dust-to-metal ra¬ 
tio. We find moderate variation in the dust-to-metal ratio 
within the galaxy, results consistent with the visualisations 
in Figures 1 and 2. 

In comparison, the models lacking growth and destruc- 
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Figure 11. Two-dimensional histogram of dust-to-gas ratio (D) versus dust-to-metal ratio {Zdust/Z) for all gas cells in the Aquarius C 
halo at 2 = 2, 1, and 0 (left, middle, and right columns, respectively) under the fiducial (FI; top row), no growth (NG; middle), and no 
destruction (ND; bottom) models. Bins are colored with the mass-weighted metallicity of constituent gas cells. In all cases, the highest 
metallicity regions tend to have high dust-to-gas ratio, but the absence of a growth mechanism for the NG model strongly suppresses 
an increase in the dust-to-metal ratio at all redshifts relative to the fiducial model. The presence of grain destruction in the NG model 
yields a peak dust-to-gas ratio more than an order of magnitude below that found in the ND model lacking a destruction mechanism. 


tion mechanisms see less diversity among gas cells. The 
no growth model has Zdust/Z < 0.5 for essentially ev¬ 
ery cell, with the typical dust-to-metal ratio decreasing to¬ 
wards z = 0. Conversely, the no destruction model sees 
-^dust /Z > 0.4 for nearly all gas cells at z = 1 and z = 0. Ad¬ 
ditionally, by z = 0, the model without destruction witnesses 
a peak dust-to-gas ratio roughly 2 dex above that obtained 
in the no growth model. These results suggest that the pres¬ 
ence of both growth and destruction mechanisms leads to 
more variation in the simulated dust content of galaxies. 

4.4 Variations in the Initial Conditions 

In the remainder of this work, we use the fiducial feedback 
and dust models to investigate the evolution of dust within 
our sample of eight Aquarius galaxies. Figure 12 shows the 
dust mass, stellar mass, SFR, and gas-phase metallicity for 
each of these galaxies as a function of redshift. These quan¬ 


tities exclude contributions from cells in the hot halo. The 
variation in galactic dust mass is largest at low redshift, 
with galaxies D and H producing dust masses of 8 x 10 8 M© 
and 3 x 10 7 M 0 , respectively, at z = 0. Changes in the dust 
content largely trace changes in the overall metal mass. We 
note that two of the galaxies with the largest dust and metal 
masses at high redshift, A and C, are known to have progen¬ 
itors that assemble their mass more rapidly than the others 
(Wang et al. 2011). The SFRs for these galaxies cover the 
range 0.4 — 7 M© yr -1 at z = 0, and the galaxy that experi¬ 
ences the largest drop in SFR from z = 1 to z = 0, galaxy E, 
also sees the largest decline in gas-phase metallicity. While 
decreased star formation produces fewer dust grains in the 
ISM, it also weakens SN-based grain destruction, ultimately 
enabling gas-phase depletion. 

To further contrast these haloes, in Figure 13 we display 
face-on and edge-on dust surface density projections for the 
Aquarius sample at z = 0. These haloes are at various stages 
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Figure 12. Clockwise from top left, a comparison of dust mass, stellar mass, gas-phase metallicity, and star formation rate for the 
Aquarius sample of galaxies with standard dust model parameters. These quantities include contributions only from the dense galactic 
disc and not the hot halo. 


of galactic disc formation. In some cases, as for halo C, a qui¬ 
escent disc of dust has formed with an abrupt drop in dust 
surface density when moving to the CGM. Other haloes, in¬ 
cluding D, F, and G, display much more spatial variation, 
especially off of the disc plane. In particular, halo G is per¬ 
turbed by a satellite galaxy at z = 0 (Marinacci, Pakmor & 
Springel 2014). A dusty disc of comparable size was found 
in recent smoothed-particle hydrodynamical simulations of 
a 10 12 Mq halo (Bekki 2015a), resembling our results for 
halo C. Furthermore, the diversity of dust surface densities 
that we see could impact the creation of synthetic galaxy 
images and spectra using spectral energy distribution mod¬ 
elling (Silva et al. 1998; Jonsson 2006; Bernyk et al. 2014; 
Snyder et al. 2015; Torrey et al. 2015). 

Figure 14 shows the projected dust-to-metal ratio of 
these haloes at z = 0 and indicates that the dust-to-metal 
ratio exhibits fluctuations within the galactic disc. Ana¬ 
lytical arguments suggest that dust-to-metal gradients can 
be used to estimate the relative strength of dust growth 
and destruction in the ISM, with flat dust-to-metal pro¬ 


files expected in the absence of dust growth and destruc¬ 
tion (Mattsson, Andersen & Munkhammar 2012). Rather 
compact galaxies, as in haloes E and H, see less spatial vari¬ 
ation in the dust-to-metal ratio than the larger galactic discs 
in Figure 13, which tend to have more pockets of low and 
high dust-to-metal ratio. Stellar density projections for the 
Aquarius suite have previously been made using the same 
physical model without dust (Marinacci, Pakmor & Springel 
2014) and show that depressions in dust-to-metal ratio are 
similar in size to regions of high stellar density. In particu¬ 
lar, haloes A, C, F, and G have some of the largest stellar 
discs and also some of the most variation in dust-to-metal 
ratio. While our dust model does not properly account for 
thermal sputtering and grain-grain collisions in hot regions 
of the CGM, Figure 14 suggests the presence of noticeable 
dust-to-metal ratio fluctuations in Milky Way-mass galactic 
discs that trace stellar structure. 

To investigate the dust-to-metal ratio evolution more 
quantitatively, in Figure 15 we plot the dust-to-gas and 
dust-to-metal ratios as a function of redshift for these galax- 
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Figure 13. Dust surface density maps for all eight simulated Aquarius haloes at z = 0 using the fiducial model. For each halo, face-on 
(top) and edge-on (bottom) projections are shown. The scale bar in the upper left indicates 25kpc. Projections were performed in a 
cube of side length 150 1 kpc centered on the halo’s potential minimum. All haloes form galactic discs in which dust is concentrated, 

although the discs vary in size. 


ies. The dust-to-gas ratios tend to evolve in similar fash¬ 
ion, reaching Milky Way-like values of (0.8 — 1.8) X 10 -2 by 
z = 0. At z = 7.5, where A1689-zDl has been observed to 
have a near-Galactic dust-to-gas ratio (Watson et al. 2015), 
the Aquarius suite of galaxies sees dust-to-gas ratios slightly 
under 10 -4 . Our results suggest that A1689-zDl and other 
dusty high-redshift galaxies do not arise from Milky Way- 
like progenitors and motivate further study of the dust-to- 
gas ratio over a range of halo masses. Recent models predict 
that larger haloes see more dramatic growth in the dust- 
to-gas ratio towards low redshift and that smaller systems 
have less dynamic dust-to-gas ratios (Bekki 2015a). A larger 
statistical sample of galaxies would also enable comparison 
to observed dust mass functions (Dunne et al. 2000, 2011). 

The dust-to-met al ratios in Figure 15 display more di¬ 
versity, with z = 0 values ranging from roughly 0.6 to 
0.7. While high depletion has been observed for some el¬ 
ements, including Si, Mg, and Fe, in the local interstellar 
cloud (Kimura, Mann & Jessberger 2003), this wouldn’t 


greatly change the overall dust-to-metal ratio. Estimates of 
the dust-to-metal ratio for typical galaxies, including the 
Milky Way and Magellanic Clouds, he closer to 0.5 (Aguirre 
1999; De Cia et al. 2013). Ah show an increase in dust-to- 
metal ratio from z = 2 to z = 0, expected of galaxies (Inoue 
2003). Before z ~ 3, about half of galaxies display a hat 
dust-to-metal ratio near 0.1. However, there is some varia¬ 
tion in when the Aquarius galaxies see the most growth in 
dust-to-metal ratio. Galaxies A and C, which assemble the 
majority of their mass more quickly than the others (Wang 
et al. 2011), show some of the earliest growth in dust-to- 
metal ratio, suggesting that accretion history may influence 
dust-to-metal ratio evolution. 

In Figure 16, we construct temperature-density phase 
diagrams for ah gas cells within halo C at z = 2, 1, and 0. We 
analyze gas cells in three different ways: by total dust mass, 
by mass-weighted dust-to-gas ratio, and by mass-weighted 
dust-to-metal ratio. We note that from z = 2 to z — 0 dust 
mass becomes increasingly concentrated in hot, diffuse halo 
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Figure 14. Projected dust-to-metal ratio in the inner disc region for the haloes shown in Figure 13. Face-on (top) and edge-on (bottom) 
projections are displayed for each halo column. The projection volume was a cube of side length 50 h -1 kpc. Figure 2 displayed higher 
resolution images for the Aquarius C halo at z = 2, 1, and 0. 


gas with T ~ 10 6 K and dense, star-forming ISM gas with 
T ~ 10 4 K near the equation of state transition. By contrast, 
the phase diagram region with T ~ 10 5 K has a total dust 
mass largely unchanged with redshift. While SDSS data do 
suggest that the dust masses of galactic haloes are similar to 
those found in galactic discs (Menard et al. 2010), the results 
in Figure 16 also motivate better modelling of dust evolu¬ 
tion in the diffuse halo where grains may travel at higher 
velocities, have lower sticking efficiencies when impacting 
gas-phase metals, and undergo grain-grain collisions. Thus, 
it is likely that the halo dust mass is overstated, especially 
at z — 0, but we still expect the presence of dust beyond the 
galactic disc. 

Dust-to-gas ratios tend to increase from z = 2 to 
z = 0 most rapidly in high-density regions. Gas density cor¬ 
relates more strongly with average dust-to-gas ratio than 
does temperature. We see a similar trend in dust-to-metal 
ratio, as dense gas above the equation of state transition 
sees its dust-to-metal ratio increase from roughly 0.2 to 0.6 
over this redshift range. In contrast, low-density gas with 


p ~ 10 -29 gem -3 barely changes in dust-to-metal ratio. This 
result is physically justified, since the dust growth timescale 
is shortest in high-density gas, and agrees with observations 
of elemental depletions increasing with gas density (Jenkins 
2009). 

In Figure 17, we display the radial dust surface density 
profiles for each of the Aquarius haloes at z — 0. There is 
fairly little evolution from z = 1 to z — 0, although the 
growth in virial radius causes Edust ~ 10 -3 M 0 pc -2 out to 
r ~ 200 kpc by z — 0. These Milky Way-sized haloes have 
dust profiles fairly consistent with observations of M31, with 
higher normalisations in the inner disc. The Aquarius haloes 
tend to show rather sharp drops in dust surface density when 
moving outside of the galactic disc: for example, halo C wit¬ 
nesses a decline of almost an order of magnitude in Edust 
near r « 20 kpc, although it is not as sharp as the drop seen 
for M31. The Aquarius haloes have dust profiles that decay 
from r ~ 20 kpc to r ~ 100 kpc in rough agreement with the 
observed scaling Edust oc r -0 ' 8 seen in SDSS data (Menard 
et al. 2010). 
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Figure 15. Comparison of dust-to-gas ratio (D) and dust-to-metal ratio (Z^ust/Z) for a variety of Aquarius galaxies with standard dust 
model parameters. The increase in dust-to-gas ratio can be as great as an order of magnitude from z = 2 to z = 0. There is significant 
variation in dust-to-metal ratio among the haloes at high redshift, though all settle near Z^ ust /Z « 0.65 at z = 0. 


In Figure 18 we compare the cumulative dust mass dis¬ 
tributions for the Aquarius galaxies with the observed distri¬ 
bution for M31 (Draine et al. 2014), normalising to the total 
dust mass within a radial distance of 25kpc from galactic 
center. We do not necessarily expect the Aquarius suite to 
have cumulative dust mass distributions identical to that of 
M31 but use these observational data for guiding purposes. 
The Aquarius haloes witness a slightly steeper enclosed dust 
mass distribution for r < lOkpc than is observed for M31, 
but the simulated and observed distributions are in rough 
agreement. In particular, we find in the Aquarius sample 
that the dust masses enclosed within lOkpc are between 
9 and 54 per cent of the total dust mass within 25 kpc. For 
comparison, the observed value for M31 is about 25 per cent 
(Draine et al. 2014) and recent smoothed-particle hydrody- 
namical simulations using a Milky Way-sized halo with dust 
growth mechanism similar to ours found that roughly 18 per 
cent of the dust mass inside of 25 kpc is contained within the 
innermost lOkpc (Bekki 2015a). 


4.5 Scaling Relations 

To better understand how realistically our sample of eight 
Aquarius galaxies produce dust, we compare against a num¬ 
ber of empirical scaling relations seen in recent observa¬ 
tions (Draine et al. 2007; Galametz et al. 2011; Kennicutt 
et al. 2011; Corbelli et al. 2012; Cortese et al. 2012; Remy- 
Ruyer et al. 2014). These scalings enable valuable predic¬ 
tions about the dust content in very metal-poor galaxies 
(Fisher et al. 2014) and high-redshift, dusty quasars (Fan 
et al. 2003; Bertoldi et al. 2003; Valiante et al. 2009). Also, 
these relations help determine the quantities that most effec¬ 


tively trace the presence of dust in galaxies. One of the most 
prominent of these scalings is the dust-metallicity relation, 
which finds a dust-to-gas ratio that scales fairly linearly with 
metallicity (Draine et al. 2007; Galametz et al. 2011; Kenni¬ 
cutt et al. 2011; Remy-Ruyer et al. 2014) and is often used 
to estimate the dust mass within a galaxy. It is therefore in¬ 
teresting to consider such scaling relations and to investigate 
any variations in our Aquarius sample. In the observational 
comparisons below, we compute ah quantities using dense 
ISM gas, as determined by the cut in temperature-density 
phase space of ah gas cells belonging to the main halo group 
given by Equation (10). This cut isolates the galactic disc 
and thus ignores the effect of dust growth in the CGM that 
is possibly too strong. 

We first show the dust-metallicity relation for the 
Aquarius galaxies in Figure 19, plotting dust-to-gas ratio 
versus gas-phase metallicity at z = 2, 1, and 0, with the 
gas-phase metallicity calculated from oxygen and hydrogen 
abundances. As in Figure 6, we compare against observa¬ 
tional data from Kennicutt et al. (2011) and Galametz et al. 
(2011). Previous models have been fairly successful at repro¬ 
ducing the dust-metallicity relation (Dwek 1998; Lisenfeld 
& Ferrara 1998; Calura, Pipino & Matteucci 2008; Bekki 
2015a; Feldmann 2015) and have investigated the nature and 
scatter of this relation at low metallicity (Remy-Ruyer et al. 
2014). In comparison to several models that see a strictly 
monotonic increase in dust-to-gas ratio and gas-phase metal- 
licity towards low redshift (Lisenfeld & Ferrara 1998; Bekki 
2015a; Feldmann 2015), the Aquarius haloes display more 
diverse redshift evolution while remaining in agreement with 
the observed dust-metallicity relation at z = 2 and z = 1. 
While the data for z = 0 capture the positive correlation be- 
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Figure 16. Phase diagram of gas temperature versus density for all gas cells in the Aquarius C halo at z = 2, 1, and 0 (left, middle, and 
right columns, respectively) where bins are colored according to the total dust mass of constituent gas cells (top row), mass-weighted 
dust-to-gas ratio (middle), and mass-weighted dust-to-metal ratio (bottom). Each row adopts its own color scale, and fiducial dust model 
parameters were used. The high-density region where gas cools is governed by an equation of state. At z = 0, the dust-to-gas and 
dust-to-metal ratios are smallest in low density gas. 


tween dust-to-gas ratio and gas-phase metallicity, they tend 
to lie above the observed relation, as in Figure 6. 

While some haloes, like A and B, display small drops in 
gas-phase metallicity from z — 1 to z — 0, others show much 
more pronounced behaviour. In particular, haloes C and E 
have some of the largest declines in SFR from z — 1 to z — 0 
in our sample, weakening grain destruction and enhancing 
depletion of gas-phase metals. Given that the temperature- 
density cut used to isolate dense gas in the galactic disc 
will minimise the influence of dust in the CGM, Figure 19 
suggests that the dust-metallicity relation can evolve in a 
non-monotonic fashion even near the galactic center. How¬ 
ever, the behaviour at z — 0 together with fairly high dust- 
to-metal ratios in Figure 15 indicate that gas-phase metals 
may be too heavily depleted at low redshift. Additionally, 
while there is evidence for a dust-metallicity relation that 
is less steep for metallicities below 8.0 (Remy-Ruyer et al. 
2014), suggesting that the global dust-metallicity trend can¬ 
not be fit by a single power law, the Aquarius galactic disc 


gas-phase metallicities lie above this region. Future simula¬ 
tions of full cosmological volumes will allow us to explore 
the dust-metallicity scaling at low metallicity. 

In Figure 20, we compare the dust mass-gas mass scal¬ 
ing seen for our Aquarius suite with observations of metal- 
rich galaxies from the Herschel Virgo Cluster Survey (Cor- 
belli et al. 2012). These observational data found the best fit 
scaling M gas oc , suggesting that the dust-to-gas ratio 

increases weakly with gas mass. That is, the presence of ad¬ 
ditional gas leads to increased star formation, which in turn 
produces enough dust to cause the dust-to-gas ratio to rise. 
We find that the Aquarius data obey a slightly steeper scal¬ 
ing of gas mass with dust mass, consistent with the results 
in Figure 15 showing fairly little variation in the dust-to-gas 
ratio in our sample. However, the Aquarius haloes tend to be 
more gas- and dust-rich than those from the Herschel Virgo 
Cluster Survey, and future work should investigate whether 
our observed scaling holds for lower mass systems. 

The Herschel Virgo Cluster Survey also investigated the 
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Figure 17. Radial profiles of dust surface density (Sdust) in the 
disc plane for all simulated haloes, plotted at z = 0. The radial 
profiles extend out to the virial radius of each halo as defined in 
Table 1. We compare with the observed dust profile for M31 (gray 
squares; Draine et al. 2014) and an M31 profile scaled by a factor 
of two (gray triangles), given that the Galactic dust content may 
vary from that of M31. The gray dashed line shows the E^ust oc 
r -o.8 sca ii n g observed in SDSS data (Menard et al. 2010), with 
normalisation adjusted to the Aquarius data. 



Figure 18. Enclosed dust mass as a function of radial distance for 
all haloes at z = 0, normalised to M(r' < 25 kpc), the mass within 
25kpc of the disc center. The normalised cumulative dust mass 
distribution observed for M31 is given by the gray dashed line 
(Draine et al. 2014). While there is some scatter in the Aquarius 
haloes, the enclosed dust mass profiles are similar to that of M31. 


Figure 19. Dust-to-gas ratio (. D) versus gas-phase metallicity 
for the Aquarius haloes plotted at z = 2, 1, and 0, with smaller 
circles denoting higher redshift. The same observational data of 
local galaxies from Kennicutt et al. (2011) and Galametz et al. 
(2011) as in Figure 6 are shown in gray. These galaxies all show 
decreases in gas-phase metallicity from z = 1 to z = 0 but fall 
close to the observed scatter at z = 2 and z — 1. 


relationship between dust-to-gas ratio and stellar mass, ob¬ 
serving the weak scaling D oc M*' 26 . This scaling is possi¬ 
bly weak because increased star formation injects more dust 
into the ISM through stellar production but also leads to 
SNe destroying dust more rapidly. In Figure 21, we compare 
the D — M* data from the Aquarius suite with the Herschel 
scaling and find the Aquarius data to lie close to the ob¬ 
served best fit. While the Aquarius haloes show a scatter in 
dust-to-gas ratio that roughly matches that seen in the Her¬ 
schel data, our haloes do not cover the full range of stellar 
masses observed by Herschel. 

In Figure 22, we compare the observed scaling between 
dust-to-stellar mass ratio and gas fraction (M gas /(M* + 
M gas )) seen in Herschel Reference Survey data (Cortese 
et al. 2012) with results from our Aquarius sample. The 
dust-to-stellar mass ratio is expected to anti-correlate with 
stellar mass and decrease from late- to early-type galaxies 
(Cortese et al. 2012). Late-type galaxies with high gas frac¬ 
tion see more efficient stellar injection of dust, while early- 
type galaxies are more influenced by grain destruction. The 
Aquarius suite is in good agreement with the observed dust- 
to-stellar mass ratio scaling, covering a range of gas fractions 
from less than 0.1 to 0.5. 

While the Aquarius haloes compare favorably to these 
dust scaling relations at z = 0, in the future we would like 
to perform similar analyses for multiple epochs. There are 
estimates of the dust mass function for z < 0.5 (Dunne et al. 
2011) and for high-redshift submillimetre sources (Dunne, 
Eales & Edmunds 2003). The relation between dust mass, 
SFR, and stellar mass has also been studied out to z = 2.5 
(Santini et al. 2014). These scalings would be interesting in 
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Figure 20. Scaling relation between dust mass (M dust ) and gas 
mass (M gas ) for the eight Aquarius haloes, with observational 
data from the Herschel Virgo Cluster Survey (Corbelli et al. 
2012). The dashed line represents the best fit scaling for the 
sample of 35 metal-rich galaxies in Virgo at z = 0.003, with 
Mg as oc M® u ^. Uncertainties in the regression parameters yield 
scaling fits given by the dotted lines. The Aquarius haloes are 
more massive than those in this observational sample and display 
a slightly steeper gas-dust scaling. 


Figure 22. Relation between dust-to-stellar mass ratio 
(Md U st/M*) and gas fraction (M gas /(M* + M gas )) for all simu¬ 
lated Aquarius haloes, shown as colored circles, with observational 
data from the Herschel Reference Survey at z ^ 0.006 provided 
as gray squares (Cortese et al. 2012). The observational data in¬ 
clude both galaxies classified as Hl-normal and Hi-deficient and 
omit upper limits. Both simulation and observation predict that 
galaxies with larger gas fraction have higher dust-to-stellar mass 
ratio. 


- 1.0 


-1.5 - 


- 2.0 - 


-2.5 - 


-3.0 


o 

o 


-1- 

A5:FI 

B5:FI 

C5:FI 

D5:FI 


—i-r 

O E5:FI 
O F5:FI 
O G5:FI 
O H5:FI 


Jr*: 


Corbelli+ (2012) 

i _ i 


9.5 


10.0 


10.5 11.0 

log(M*/M 0 ) 


11.5 12.0 


Figure 21. Dust-to-gas ratio (D) plotted as a function of stellar 
mass (M*) for our sample of Aquarius haloes, with the same 
Herschel Virgo Cluster Survey data as in Figure 20 given in gray 
(Corbelli et al. 2012). The observed weak scaling with D oc M* -26 
is indicated by the dashed line. The scatter of Aquarius haloes is 
similar to that in these observations. 


part because our Aquarius haloes displayed more diversity 
at high redshift than at z — 0. 

Especially at high redshift, it is important to under¬ 
stand the relative dust production strengths of SNe and 
AGB stars, since this helps determine whether dusty galax¬ 
ies and quasars could have dust content driven by stellar 
sources as opposed to interstellar dust growth. Chemical 
evolution models suggest that AGB stars played an impor¬ 
tant role in forming the large dust mass of SDSS J1148+5251 
at z = 6.4, possibly contributing over half of the observed 
dust (Valiante et al. 2009). Observations of dust in the ejecta 
of SNe fall short of the 1 M 0 or more a SN needs to produce 
in order to explain the dust content of SDSS J1148+5251 
(Todini Ferrara 2001; Sugerman et al. 2006; Dwek, Gal¬ 
liano <V Jones 2007; Lau et al. 2015), and if SNe are not 
dominant producers of dust, high SFRs might be needed to 
form dusty z > 5 galaxies (Morgan <V Edmunds 2003; San- 
tini et al. 2014). Given that the progenitors of AGB stars 
take longer to evolve off of the main sequence than SNe, it is 
unclear whether AGB stars can have a large impact at high 
redshift. In Figure 23, we compare the dust mass contribu¬ 
tions of SNe la, SNe II, and AGB stars for the Aquarius C 
galaxy as a function of redshift. We see that 80 per cent or 
more of dust originated in SNe II for z > 5, with the contri¬ 
bution from AGB stars peaking at roughly 40 per cent for 
z ~ 2. Recent investigation of the Small Magellanic Cloud 
suggests that even at low redshift AGB stars are not domi¬ 
nant producers of dust (Boyer et al. 2012). In this Aquarius 
galaxy, we see that by z = 0 about 20 per cent of dust has 
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Figure 23. Evolution of dust mass contributions from different 
stellar types for the Aquarius C galaxy using the fiducial dust 
model. The bottom subpanel shows the fraction of total dust mass 
(/tot) contributed by each component. As noted in Section 2.2.3, 
the relative proportions of dust produced by AGB stars and SNe 
are kept constant during changes in dust mass from dust growth 
or destruction. Type II SNe are responsible for roughly 80 per 
cent of dust formed in the first Gyr and two-thirds of dust at 
z = 0. This figure captures the delayed production of dust by 
AGB stars owing to their long lifetimes. 

its origins in AGB stars. While future work is needed to in¬ 
vestigate how these stellar contributions vary with galaxy 
stellar mass and metallicity, in our sample of Milky Way¬ 
sized galaxies we find SNe to be the dominant producers of 
dust. 


5 DISCUSSION AND CONCLUSIONS 

We have implemented a first very basic dust model in the 
moving-mesh code AREPO, adding to the existing galaxy for¬ 
mation physics. Our dust model accounts for local stellar 
production of dust, growth in the ISM due to accretion of 
gas-phase metals, destruction via SN shocks, and dust driven 
by stellar feedback winds. Dust is also passively advected be¬ 
tween gas cells. We track dust in individual chemical species 
as well as follow contributions from AGB stars, SNe la, and 
SNe II. 

Using this dust model, we performed cosmological 
zoom-in simulations of the Aquarius suite of eight haloes 
with the goal of understanding how dust forms in Milky 
Way-sized systems. After investigating the effect of feedback 
and the strength of various dust model components, we used 
the full sample of Aquarius haloes to compare to a range of 
dust observations at low redshift. We summarise our main 
findings as follows. 

(i) Variations in stellar and AGN feedback impact the 
predictions made by our dust model. The absence of any 


feedback led to a z — 0 dust-to-metal ratio of 0.2 for one 
Aquarius galaxy, roughly a factor of three lower than the 
value predicted by our fiducial feedback model. Exclud¬ 
ing feedback led to high star formation and more efficient 
supernova-based dust destruction. On the other hand, a 
model with fast stellar feedback-driven winds helped sweep 
dust away from the galactic center and enhance depletion 
of gas-phase metals due to weakened star formation, raising 
the halo-wide dust-to-metal ratio to over 0.7. 

(ii) Allowing gas-phase metals to deplete onto dust grains 
and not contribute to metal-line cooling has the potential 
to affect gas-phase metallicity. We see a difference in halo 
gas-phase metallicity between runs with and without dust 
of roughly 0.5 dex at z = 0. This result also motivates the 
future inclusion of dust cooling channels. 

(iii) In the absence of dust growth, dust-to-gas and dust- 
to-metal ratios are suppressed by more than an order of 
magnitude at z — 0 from the values predicted by our fidu¬ 
cial model with dust growth and destruction mechanisms. 
Without growth, we obtain a dust-to-gas ratio of 3 x 10 -4 
for one Milky Way-like galaxy, below the Galactic value of 
roughly 1C) -2 , and a dust-to-metal ratio Zdust/^ < 0.05. The 
differences between the fiducial model and no dust growth 
model are noticeable for z < 5. 

(iv) For z > 5, the dust mass and dust-to-gas ratio 
produced by a model adopting a constant dust destruc¬ 
tion timescale of ra = 0.5 Gyr are largely indistinguish¬ 
able from those produced by a model with no dust destruc¬ 
tion. By z = 0, the constant destruction timescale model 
yields results similar to those from the fiducial model with 
local SNe-based dust destruction. At z = 0, the model 
with constant destruction timescale produces a dust mass 
AGust = 3 x 10 8 M© in the galactic disc with a dust-to-gas 
ratio D = 10 -2 . 

(v) The Aquarius haloes form dusty galactic discs with 
typical surface densities Edust ~ 10 -1 M©pc -2 . The inclu¬ 
sion of thermal sputtering and grain-grain collisions would 
likely reduce the dust-to-metal ratios seen in the CGM. 
Dust-to-gas ratios for these Milky Way-like galaxies increase 
by roughly an order of magnitude from z = 3 to z = 0 and 
cover the range (0.8 — 1.8) x 10 -2 at z = 0. 

(vi) The predicted ISM dust content of the Aquarius 
haloes is consistent with a number of observed scaling re¬ 
lations at z = 0, including scalings between dust mass 
and gas mass, dust-to-gas ratio and stellar mass, and dust- 
to-stellar mass ratio and gas fraction. While the overall 
trend of the dust-metallicity relation is reproduced over 
8.1 < 12 + log(0/H) < 9.0 at z = 2 and z = 1, the Aquarius 
galaxies do not evolve along strictly monotonic tracks from 
z = 2 to z = 0. 

(vii) Our dust model tracks contributions from individual 
chemical species, but additional physics is needed to cap¬ 
ture differences between silicate and graphite grain types. As 
shown in Figure 24, currently silicate and graphite grains are 
distributed in essentially the same manner. While we may 
expect slight variations, given that AGB stars and SNe pro¬ 
duce dust of different compositions and the dust originating 
from SNe may be more readily destroyed in SN shocks, fu¬ 
ture work should include starlight-driven radiation pressure 
on dust grains. Given that the dust composition will im¬ 
pact the grain size distribution (Mathis, Rumpl V Nordsieck 
1977; Kim, Martin V Hendry 1994) and radiation pressure 
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Figure 24. Surface densities of silicate (top) and graphite (bot¬ 
tom) components of dust for the Aquarius C halo at z = 2, 1, 
and 0 (left, middle, and right columns, respectively). The details 
of the projections are the same as those in Figure 1. The sil¬ 
icate and graphite dust components largely trace one another, 
although future work accounting for radiation pressure and grain 
size distributions may yield more diversity. 

is suggested to be important in simulations treating dust 
(Murray, Quataert V Thompson 2005; Bekki 2015b), this 
physics may help differentiate between silicate and graphite 
grains. 

While our current dust model is fairly passive and sim¬ 
plistic, it reproduces a number of observed dust scaling rela¬ 
tions surprisingly well at low redshift and is a step towards 
a more complete treatment of dust in cosmological galaxy 
formation simulations. In the future, we can explore a live 
dust model that includes dust-gas interactions more directly, 
enables dust cooling channels, and allows for radiative forces 
to affect the motion of grains. Together with simulations of 
full cosmological volumes, this work will lead to a better un¬ 
derstanding of the full diversity of dust seen in the Universe. 
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